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Abstract. A new class of piecewise linear methods for the numerical solution of 
the one-dimensional Euler equations of gas dynamics is presented. These methods are 
uniformly second-order accurate, and can be considered as extensions of Godunov’s 
scheme. With an appropriate definition of monotonicity preservation for the case of 
linear convection, it can be shown that they preserve monotonicity. Similar to Van 
Leer’s MUSCL scheme, they consist of two key steps: a reconstruction step followed by 
an upwind step. For the reconstruction step, a monotonicity constraint that preserves 
uniform second-order accuracy is introduced. Computational efficiency is enhanced by 
devising a criterion that detects the ‘smooth’ part of the data where the constraint is 
redundant. The concept and coding of the constraint are simplified by the use of the 
median function. A slope-steepening technique, which has no effect at smooth regions 
and can resolve a contact discontinuity in four cells, is described. As for the upwind 
step, existing and new methods are applied in a manner slightly different from those 
in the literature. These methods are derived by approximating the Euler equations via 
linearization and diagonalization. At a ‘smooth’ interface, Harten, Lax, and Van Leer s 
one-intermediate- state model is employed. A modification for this model that can 
resolve contact discontinuities is presented. Near a discontinuity, either this modified 
model or a more accurate one, namely, Roe’s flux-difference splitting, is used. The 
current presentation of Roe’s method, via the conceptually simple flux-vector splitting, 
not only establishes a connection between the two splittings, but also leads to an 
admissibility correction with no conditional statement, and an efficient approximation 
to Osher’s approximate Riemann solver. These reconstruction and upwind steps result 
in schemes that are uniformly second-order accurate and economical at smooth regions, 
and yield high resolution at discontinuities. 
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1. Introduction. Among methods for the numerical solution of conservation laws 
in general, and the Euler equations in particular, upwind schemes are very popular due 
to their accurate shock- capturing capability. Review articles by the pioneers of these 
methods are readily available (Harten, Lax, and Van Leer 1983, Roe 1986, 1989). Sev- 
eral formulations of upwind schemes, which generally extend that of Godunov (1959), 
exist [see also Boris and Book (1973)]. Here, we essentially employ the projection- 
evolution approach advocated by Van Leer (1979, 1985). In this approach, an upwind 
scheme consists of two key steps: a reconstruction (projection) step where, within each 
cell, the data axe approximated by polynomials; and an upwind (evolution or approxi- 
mate Riemann solver) step where the average fluxes at each interface are evaluated by 
a procedure that takes into account the ‘wind’ (or ‘wave’) directions. The Godunov 
scheme employs the simplest reconstruction — the piecewise constant function — and 
the most accurate upwind procedure — the exact solution of the Riemann problem. Al- 
though this method is only first-order accurate and computationally costly, it is the 
basis of most modem upwind schemes. A great deal of effort has been spent to enhance 
•accuracy by using higher-order polynomials in the reconstruction step, and to improve 
efficiency by approximating the solution of the Riemann problem in the upwind step. 

First, we discuss the reconstruction step. For simplicity, we consider only linear 
reconstructions here. Accurate interpolations (linear and higher order) are derived by 
assuming that the data are smooth. If the data have a discontinuity, these interpo- 
lations always create oscillations in the solutions. To prevent such oscillations, Van 
Leer (1974, 1977, 1979) introduced a monotonicity constraint for his MUSCL scheme: 
the reconstruction must preserve the monotonicity of the data. At a smooth region 
with nonzero slope, the constraint is redundant provided the mesh is fine enough, and 
the MUSCL scheme is second-order accurate. At a discontinuity, the constraint takes 
effect to prevent oscillations. A major drawback of the constraint is that it also takes 
effect adjacent to an extremum and causes accuracy to degenerate to first order. In 
other words, high resolution at discontinuities is obtained at the expense of accuracy 
near extrema. The popular total variation diminishing (TVD) schemes (Harten 1983, 
Sweby 1984, Roe 1985) share the same advantages as well as disadvantages since they 
satisfy Van Leer’s monotonicity constraint. 
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Harten and Osher (1987) solved the loss of accuracy problem by introducing the 
nonoscillatory piecewise linear reconstruction, which involves no constraint. Their 
UNO scheme is unif ormly second-order accurate and does not create oscillations, but 
it is diffusive, e.g., it smears contact discontinuities. The reason for this smearing effect 
is that, roughly speaking, the slope of the UNO reconstruction is closest to zero among 
all the slopes that are second-order accurate {0(h 2 )); therefore, its error is one of the 
largest in this class (Huy nh 1993). Modifications of the UNO and ENO (Harten et 
al. 1987) schemes to improve the resolution at contact discontinuities can be found 
in Harten (1989), Yang (1990), and Mao (1991). These modifications, however, are 
somewhat complicated. For more developments of these schemes, see Shu and Osher 
(1989), Rogerson and Meiburg (1990), and Shu (1990). For a technique that senses an 
extremum and turns off the constraint there, see Leonard and Niknafs (1991). 

In Huynh (1989, 1993), the author presented a geometric framework in which Van 
Leer’s constraint can be combined with Harten and Osher’s nonoscillatory interpola- 
tion. For scalar conservation laws, the resulting SONIC schemes are uniformly second- 
order accurate and nonoscillatory. In addition, each TVD scheme corresponds to a 
SONIC counterpart, and the UNO scheme is the most diffusive member of the SONIC 
class. 

In this paper, the author’s approach is extended to the system of Euler equations. 
Similar to the UNO scheme, we employ a five-point stencil for the reconstruction. 
While the UNO slope is ‘closest’ to zero, we start with the most accurate one for this 
stencil: the slope of the quartic through these five points. Next, since monotonicity 
constraints are derived for the convection equation, and characteristic variables are 
the quantities that are being convected, the constraints are applied to these variables. 
Our constraint consists of two limits (bounds): the lower one, which assures uniform 
second-order accuracy, is defined by a slope similar to that of the UNO scheme, the 
upper one, which controls oscillations, is derived by making use of the upper limit 
for the MUSCL scheme. The constraint requires the final slope to lie between these 
two limits. This requirement is conveniently enforced by using the median function: 
the final slope is the median of the quartic slope and the above two limits. As such, 
the algorithm is costlier than that of the UNO scheme. To save computing time, we 
present a simple criterion that detects the ‘smooth’ part of the data: if a cell is in the 
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‘smooth’ region, then the monotonicity constraint has no effect, and the slope reduces 
to the quartic formula for the primitive variables. With this criterion, the characteristic 
variables and the constraints are necessary only for about six to eight cells around each 
discontinuity. The resulting algorithm for the reconstruction step is therefore accurate 
and efficient; however, it still smears a contact discontinuity, although to a lesser extent 
them the MUSCL or UNO schemes. To improve the resolution in this situation, we 
present a simple slope-steepening technique that has no effect at the smooth part of 
the data and generally can resolve a contact discontinuity in four cells or less. 

Once the slopes (space partial derivatives) are defined, the time partial derivatives 
can be calculated by the non-conservation form of the Euler equations (Harten et 
al. 1987, Roe 1989). At each interface, we can march half a time step via Taylor series 
expansions, but the two adjacent cells yield two sets of results (states); we need to derive 
one set of fluxes from these two states. This is accomplished by one of the upwind 
procedures discussed below. Note that this half-step calculation is a simplification 
of an observation by Hancock (Van Albada, Van Leer, and Roberts 1982). It helps 
simplify our algorithm at the ‘smooth’ part of the data (below); in addition, it makes 
the treatment of source terms straightforward. 

While the reconstruction step is mainly numerics, the upwind step is where physics 
plays its role. There are several upwind procedures in the literature (Roe 1981, Steger 
and Warming 1981, Osher and Solomon 1982, Van Leer 1982, Harten, Lax, and Van 
Leer 1983, Roe and Pyke 1984, Roe 1986, Einfeldt 1988, Toro 1991, Liou and Steffen 
1993). The most popular method is perhaps that of Roe, but the perfect one is yet 
to be found (Roe 1989, Roberts 1990). Upwind procedures are generally perceived 
as difficult to formulate and implement. In this paper, existing and new methods 
employed are derived from a simple and unified point of view: the Euler equations 
in conservation form are approximated by linearization and diagonalization. For the 
flux-vector splitting introduced below, what is crucial is that the fixed state for the 
linearization is chosen so that the resulting fluxes depend continuously on the data, 
and the left and right flux vectors are diagonalized by the same matrices, i.e., those 
formed by the eigenvectors of the fixed state. 

Since our reconstruction step already distinguishes between the ‘smooth’ and ‘non- 
smooth’ parts, a different upwind flux can be applied to each part. At a ‘smooth’ 
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interface, the left and right states are close to each other, and a simple model is suffi- 
cient; in fact, a sophisticated one does not significantly improve the results. Here, we 
make use of Harten, Lax, and Van Leer’s one-intermediate-state model, which obtains 
the flux via a scalar calculation. Combining with the quartic formula, at a smooth 
region, the computing effort is therefore only slightly more than that of a central- 
difference scheme with artificial viscosity (Jameson, Schmidt, and Turkel 1981). Using 
information from the reconstruction step, we also present a modification for the one- 
intermediate-state model that can resolve contact discontinuities. Near a discontinuity, 
either this modified model or a more accurate one, namely, Roe’s flux-difference split- 
ting (1981), is employed. Our presentation of Roe’s method as a flux-vector splitting 
(Steger and Warming 1981) has the advantage of conceptual simplicity. It also estab- 
lishes a connection between the two splittings. This connection, in turn, facilitates the 
derivation of an admissibility (or entropy) correction with no conditional statement 
for the flux- difference splitting. In addition, it leads to an efficient approximation to 
Osher’s approximate Riemann solver (Osher and Solomon 1982), which is named the 
simple- wave upwind flux. 

At smooth regions — and generally, most of the flow field is smooth the scheme 
is accurate and economical due to the quartic formula and the one-intermediate-state 
model. Near a discontinuity, the monotonicity constraint takes effect to prevent oscil- 
lations without causing a loss of accuracy. Moreover, the upwind step for this case, 
via either Roe’s flux-difference splitting with admissibility condition or the simple-wave 
upwind flux, yields high resolution at shocks as well as contact discontinuities, and non- 
physical solutions are also excluded. The resulting schemes are efficient and, as will be 
shown by numerical experiments, their accuracy compares favorably with higher-order 
schemes. 

This paper is essentially self-contained and is organized as follows. In §2, we review 
the various forms of the Euler equations and the corresponding characteristic variables. 
The discretization and the half time step is carried out in §3. The two main sections 
4 and 5 are independent of each other: §4 deals with the reconstruction step; §5, the 
upwind step. The algorithm is extended to the case of flows with area variations in 
§6. Numerical experiments are shown in §7. Finally, the discussion and conclusions 
axe presented in §8. 
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2. Forms of the Euier equations of gas dynamics. The Euler equations axe 
derived from the physical principles of conservation of mass, momentum, and energy. 
In this form, called the conservation form, the equations axe valid whether the flow is 
smooth or has discontinuities. To assure correct shock speed (Lax 1954), therefore, we 
solve the equations in conservation form. The problem thus reduces to a flux estimate 
at each interface. To calculate this flux, we make use of the systems of equations in non- 
conservation form for both the primitive and conservative variables. The characteristic 
vaxiables associated with each of these two systems also play an essential role: upwind 
algorithms axe generally derived for a scalar convection equation; for systems, the 
quantities that are being convected are the characteristic vaxiables. We review the 
various forms of the Euler equations in this section [see also Hirsch (1990)]. Different 
approaches that make use of these forms can be found in Roe and Pyke (1984), Roe 
(1986, 1989), and Toro (1991). 

2.1. Conservation form. The one-dimensional flow of an inviscid and compress- 
ible gas obeys the conservation laws for mass, momentum, and energy: 


dU dF(U) 
dt + dx 


= 0 , 


( 2 . 1 a) 


U = 



m 


F = I ( m 2 /p) + p 

(e+p)m/p l 


(2.16, c) 


where t is time, x distance, p density, m momentum, e total energy per unit volume, 
and p pressure. Let 7 be the ratio of specific heats, then for a perfect gas, the system 
is completed by the equation of state 


P = ( 7 - l)[e-m 2 /(2p)]. 


(2.2a) 


Denote by u the velocity, then m = pu, and the above yields 


P _1 2 

e = + -pu . 

7 — 1 2 

The flux vector is more commonly written as 

\{e + p)uj 


(2.26) 


(2.2c) 
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the above expression is also convenient for the purpose of coding. 

2.2. Primitive variables. With superscript T representing the transpose, let V 
be the column vector of primitive variables, 

V = (p,u,p) T . (2.3a) 



One can easily verify that the columns of the following matrix are the right eigenvectors 
corresponding to the above eigenvalues, 

Rp= ( -c/p 0 c/p | . (2.6a) 

V c 2 0 c 2 ) 

Now set 

L P = Rp 1 ; (2.66) 

then (o -/./(2c) 1/ (2c 2 ) \ 

Lp = 1 0 -1/c 2 , (2-6c) 

\0 p/ (2c) l/(2c 2 ) / 

and the rows of L p are the left eigenvectors of A p . It follows from the definitions of Rp 
and L p that 

LpApRp = A, A p = RpALp, (2.7a, b) 


7 



where 


/ AW 0 0 \ 

A = 0 A< 2 > 0 . (2.7c) 

V 0 0 A ( 3 >/ 

Note that we do not need the subscript p for A because, as shown below, this matrix is 
identical to that of the conservative variables. Also observe that the matrices L p and 
R* axe functions of p and c only (or p and p); they do not depend on u. 

Next we define the characteristic variables associated with V. Let V = (p,u,p) T 
be a fixed state. Let A p , L p , R,, and A be the corresponding matrices given by (2.3c), 
(2. 6c, a), and (2.7c). By linearizing (2.3b) at V, 


dV 


dV 


JT + Ar ~fc-°- 


(2.8a) 


Since A is a constant matrix, the above equation is linear; it is an approximation to 
the nonlinear equation (2.3b) in a neighborhood of V provided that the. solution V is 
smooth. By (2.7b), 

A p — RpAL P . 

In addition, because L p is a constant matrix, and it is the inverse of Rp, equation (2.8a) 
is equivalent to 

/if. V . AT . V 

( 2 . 86 ) 


^ v + A^X=o. 


dt ' dx 
We can now define the characteristic variables. Set 


W p = L p V ; 


(2.9a) 


then it follows from (2.6c) that 


W p = 


and (2.8b) takes the form 


/ -[ / »/(2c>] « + [1/(2^)] p' 
p- [l/£ ! ]p 
[^/(2c)]u+[l/(2c 2 )]p 


aw f + A aw, = 0 . 


(2.96) 


(2.9c) 


dt dx 

Because A is a diagonal matrix, the above system is decoupled, and each component 
corresponds to a convection equation of speed u — c, u, and u + c, respectively. For 
this reason, the components of W p axe called characteristic variables. Note that all 
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components of W p have units of density; as a result, a comparison among them makes 
sense. The second component corresponding to the speed u is called the linear char- 
acteristic; the other two are the nonlinear ones. The above characteristic variables are 
simpler than those defined from the conservative variables below. 


2.3. Conservative variables. At smooth regions of U, equation (2.1) can be 
written in non-conservation form via the chain rule: 


au , dU „ 

ar + Ae ar- 0 ’ 

(2.10a) 

dF (dF dF dF\ 
dU \ dp ’ dm ’ de ) 

(2.106) 


where A c is the Jacobian matrix, 

A c = 

A straightforward calculation using (2.2a) yields 

( 0 1 

A c = (7 - 3)u 2 /2 (3 - 7)u 

\ (7 — l)u 3 — 7ue/p -3(7 - l)u 2 /2 + 7e/p 

To diagon aliz e A c , we make use of the transformation from the primitive variables 
to the conservative ones (Warming, Beam, and Hyett 1975). This transformation will 
also play a key role in the relations between the flux- vector and flux-difference splittings 
later. Denote the Jacobian matrix of this transformation by M: 



(2.10c) 


M 


dU 

dV 


k dp 


dU dU \ 
’ du ’ dp ) 


(2.11a) 


Then it follows from the definitions of U and V and (2.2b) that 

/ 1 0 

M = 


0 

up 0 

\u 2 / 2 pu 1/(7 — !) 


(2.116) 


As a result, 


/ 


M" 1 = 


1 

-u/p 


0 

1 Ip 


0 

0 


(2.11c) 


\ (7 - 1 )u 2 /2 -(7 - 1)« 7 ~ 1 . 

Next, we relate A c and A p by using M. Apply the chain rule to (2.10a), 


( 2 . 12 ) 
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As a consequence of the above and (2.3b), 


A„ = M _1 A r M. 


This fact can also be verified by carrying out the right hand side matrix multiplications 
using (2.10c) and (2.11b, c); the result is identical to (2.3c). Substitute (2.13) into 
(2.7a), 

(L p M -1 )A c (MRp) = A. (2.14a) 


Now set 


Then 


L c = L P M ', R- = MRp. 


Lc.AcR'c — A, 


(2.146, c) 


(2.15) 


the matrices L c and R c are the inverse of each other; the columns of R c are the right 
eigenvectors of A c , and the rows of L c , the left. Expressions (2.14c), (2.11b), and (2.6a) 
yield 

1 1 1 \ 

R c = u — c u u + c I , (2.16) 

{H — uc u 2 / 2 H + ucJ 


where H is the total enthalpy, 


_ e + p _ u 2 7P + — 

P 2 (7-1 )p 2 7 -I' 


(2.17) 


Similarly, from (2.14b), (2.11c), and (2.6c), 

/ M/2 + (7 - l)M 2 /4 [-1 - (7 - l)M]/(2c) (7 - l)/(2c 2 ) 

_ l-(7 -l)M 2 /2 (7 — 1)M/ c -(7 -l)/c 2 

-M/2 + (7 - l)M 2 /4 [1 — (7 — l)M]/(2c) (7 — l)/(2c 2 ) 


(2.18) 


where the Mach number M is defined by 

M = u/c. (2.19) 

Later, we will make use of the fact that L c and R c depend only on any two of the three 
variables u, H , and M; they do not depend on p. In fact, the matrix A c can also be 
written in a form that only depends on u and H: by (2.17), 

p = pH- e; (2.20) 
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substituting the above into ( 2 . 2 b), 


2£ = i ,+ lzi u ’ ; (2.21) 

p 2 

as a result, 

0 1 0 \ 

A c = (7 - 3)u 2 /2 (3 - 7)« 7 ~ 1 • ( 2 - 22 ) 

^ (7 — l)u 3 /2 — uH — ( 7— 1 )u 2 + H 7 u ) 

The derivation and definition of the characteristic variables associated with U are 
similar to those for V in (2.9): with L c defined by (2.18) at a fixed state U, set 

W c = L C U; (2.23a) 


then it follows from (2.10a) and (2.15) that 


dW c , idW c 

^r + A -ar 


= 0 . 


(2.236) 


By ( 2 . 6 c), L p is sparse; by (2.18) L c is dense. Consequently, W p is employed whenever 
possible, and W c is used only when necessary. 

Before moving to the next section, we make the following remark. At smooth regions 
of the solution, the three forms (2.1), (2.3), and (2.10) axe equivalent; at a discontinuity, 
however, only the conservation form ( 2 . 1 ) is valid. In each cell, therefore, since the 
solution is assumed to be smooth, we can employ (2.3), i.e., the primitive variables, 
for the purpose of reconstruction because it is the simplest of the three. Of course the 
conservative variables ( 2 . 10 ) can be used for this step at the cost of extra computing 
time. The results for these two reconstructions are essentially identical for the second- 
order methods presented below. (Note that for higher-order methods, the advantage of 
employing ( 2 . 3 ) over ( 2 . 10 ) is no longer clear because deriving accurate averages for the 
primitive variables is more involved.) At an interface, however, there is a discontinuity; 
consequently, we employ ( 2 . 1 ) to obtain the flux. 


3. Discretization. Let Xj = jh, j = 0, 1,2, . . . be a uniform mesh, xj +1 / 2 = 
(j + 1/2)6 be the boundary between the jth and j + 1 cells, and t n = nr be the time 
level; here, h is the mesh size, and r, the time step. Let U” be an approximation to 
the average value of the solution vector U in the j- th cell at time t n , 

U!- « j f 3+ ' U(x,f n )dx, (3.1) 
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where the integration acts on each component of U. Assume that we know U* , for all 
j; we want to calculate U/ +1 , where the time step r satisfies the CFL condition 


|max 




(3.2) 


In this and the next two sections, we only consider the interior points; boundary 
conditions employed in this paper are standard and will be addressed as needed. 

To derive second-order methods, we use the midpoint rule: 

j\ n-|-l/2 jpJi+l/2 


U? +1 - U| 


+ 


• j+1/2 


^ = 0 , 


(3.3a) 


or equivalently, 


u? +1 = u* + - F£$). 


h 


(3.36) 


For simplicity of notation, we omit the superscript n when there is no confusion; e.g., 
Uj denotes U”; the superscripts n + 1/2 and n + 1, however, are retained. 

The problem, therefore, reduces to estimating from the data {U,} (the 

braces represent all j). To obtain these fluxes, we first calculate the primitive variables 
{V/} in a straightforward manner by their definitions. Next, in each jth cell, and for 
t n < t < f n+1 , we approximate V (x, t) by a linear function Rj(x, t): 


Rj(z, t) = Vj + (x - xy)S j + (t - t n )Tj, (3.4) 


where Sj and Tj respectively approximate V r (xj,f n ) and V t (xj, t n ). Suppose, for the 
moment, Sj is known. We march to time f n+1 / 2 as follows. Since Vj is also known, Tj 
can be estimated by using the equations for the primitive variables (2.3): 


Tj = — (Ap)j- Sj. (3.5) 

At time t n+1 ^ 2 and at the two interfaces of the jth cell, 

Rj(*,--i/ 2 .i“ +l/2 ) = Vj - |s,. + jT„ (3.6a) 

R,'(z i+1/2 , ("«/’) = Vj + |s, + jTj. (3.66) 

At each interface j + 1/2, we now have two values for the primitive variables: one from 
the Taylor series expansion in the jth cell, namely, Rj(xj + i/ 2 , f n+1 / 2 ); and one from 
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that in the j + 1 cell, Ry +1 (x ;+ i/ 2 , £ n+1/2 ). From these two values, we will calculate the 
flux via one of the upwind procedures that employ the conservation form (2.1). 

It is only here that the ‘wind’ (or ‘wave’) directions come into play. 

As a result of the above arguments, the problem reduces to a slope estimate and 
an upwind procedure for the flux. We remark that these two steps are independent of 
each other. 

Note that the calculation of the time partial derivative from the spatial one via (3.5) 
was employed in the derivation of the Lax-Wendroff (1964) scheme. In the context of 
reconstruction, this procedure for the conservative variables was observed by Hancock 
(Van Albada, Van Leer, and Roberts 1982). It can be traced back to the Cauchy- 
Kowalewsky theorem (Courant and Hilbert 1962, Harten et al. 1987). 

4. Slope estimates. This section deals only with interpolations and constraints; 
the ‘wind’ direction does not play any role. When the data are smooth, accurate 
formulas for the slope are readily available (parabolic or quartic). Near a discontinuity, 
however, such formulas are no longer accurate; in fact, they yield a poor representation 
of the data and create oscillations in the solutions. By imposing various constraints on 
the slope, a better representation of the data can be obtained; the resulting schemes 
are generally also free from oscillations. The constraints should be designed so that 
at smooth regions of the data, they have no effect; i.e., they do not alter the original 
slope. The monotonicity constraint presented below consists of two limits (bounds): a 
lower one, which is 'closest’ to zero among all slopes that are accurate to 0{h ), and an 
upper one, which prevents the slope from becoming too steep. The constraint requires 
that the final slope lies between these two limits. 

4.1. Constraints for the slopes. To be precise, we need a few definitions. Let 
zi, . . . , Zfc be real numbers, and let 

a = min(zi, . . . , z k ), P = max(zi,...,z fc ). (4.1) 

We define I [z u ...,z k ] to be the smallest closed interval containing z u . . . , z k ; i.e., 

I [ z i ■>•••■> Zfc] = [<*,£]■ 

Next, let the median of three numbers be the one that lies between the other two. 
Then it follows that the median lies in the interval defined by any two of the three 
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arguments; e.g., if x, y, and z are real numbers, then 

median (x, y,z) € /[r,y], median (z, y, z) el [y, z\. (4.2) 

It is very convenient to enforce a constraint by using the median function. As an 
example, suppose y\ and y 2 are the two limits for a quantity y; i.e., the constraint is 
that y lies in the interval I [yi, y 2 ]. Then this constraint is satisfied by y/ (/ for final) 
defined by 

y/ = median (y, y 1} y 2 ). (4.3) 

Note that if y lies between yi and y 2 , then y/ = y, and the constraint does not alter the 
original estimate. If y lies outside the interval /[yi,y 2 ], then between the two values 
yi and y 2 , yj equals the one that is closer to y. 

We now make an observation concerning accuracy: if two of the three arguments 
tire accurate to a certain order, then the median yields a value that is accurate to the 
same order. More precisely, suppose x and y approximate an exact value y to 0(h k ); 
then median (r, y,z) is also accurate to 0(h k ) regardless of the value z: 

median (x, y, z) = y + 0(h k ). (4.4) 

Indeed, since x and y are accurate to 0(h k ), so is any value between x and y. By (4.2), 
median (x, y, z) lies between x and y; as a result, it is accurate to the same order. 

With zi, . . . , Zfc being real numbers, let minmod (zi, . . . , z*) be the function which 
returns the smallest argument if all arguments are positive, the largest if they are all 
negative, and 0 otherwise. This definition is equivalent to 

minmod (zi, . . . , Zfc) = median (or, /?, 0) (4.5a) 

where a and 0 are given by (4.1). For the purpose of coding, 

minmod (zi,..., Zfc) = smin(|zi|, . . ., \z k \) (4.56) 

where 

s = |[sgn (zO + sgn (z 2 )] §|sgn (zi) + sgn (z 3 )| . . . ||sgn (z x ) + sgn (z fe )|, (4.5c) 

and sgn (z) = 1 if z is positive; sgn(z) = — 1 if z is negative. In the above formula 
for s, the expression within the square brackets yields the appropriate sign. Note that 
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there is only one set of square brackets, but there are several sets of absolute value 
sign. Clearly, if any two arguments are of opposite sign, then one of them must be of 
opposite sign with Z\\ consequently, s = 0. In general, s equals either —1, 0, or 1. Also 
observe that if any argument z t - equals 0, then it does not matter whether sgn (z t ) is 
defined as —1, 0, or 1 in (4.5c) since the resulting minmod, by (4.5b), is 0. 

With two arguments x and y, minmod (x, y) is simply the median of x, y, and 0. 
In this case, expressions (4.5b,c) reduce to 

minmod (x,y) = | [sgn (x) + sgn (y)] min(|x|, |y|). (4.6) 


The minmod function is more commonly defined as 


minmod (x, y) = f “f* >0 ’ 

v lo otherwise 

(e.g., Harten 1983, Sweby 1984, Haxten and Osher 1987). 

For convenience of programming, we express the median function in terms of min- 
mod: 


median (x, y, z) = x + minmod (y — x, z — x) 

= y + minmod (x — y,z — y) (4-7) 

= x + \ [sgn (y - x) + sgn (z - x)] min(|y - x|, |z - x|). 


Let the extended-minmod function xm (x, y) be defined by 


xm(x, y) = median (x,y, — x — y). (4.8a) 

(In terms of Sweby’s limiter function (1984), xm(x,y) = x median(l,r, —1 — r) where 
r = y/x.) It follows from (4.7) that 

xm (x, y) = -(x + y) + |[sgn(2x + y) + sgn(2y + x)] min(|2x + y|, |2y + x|). (4.85) 

As explained below, the above function is slightly different from the minmod func- 
tion. If x and y are of the same sign, then the two functions return the same value, 
namely, the one with smaller modulus. If x and y axe of opposite sign, however, the 
minmod function returns 0, while the extended-minmod returns the value with smaller 
modulus provided that |x/y| > 2 or |x/y| < 1/2 [see also expression (2.27) in Huynh 
(1993)]. Loosely speaking, the extended minmod function returns the argument of 
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smaller modulus nearly all the time, and next to a discontinuity, it preserves the slope 
of the smooth side more often than the minmod function. 

We axe now ready for the constraint. Since monotonicity constraints axe derived 
for the convection equation (Van Leer 1977), and the characteristic variables are the 
quantities that axe being convected, we interpolate W p , the characteristic variables 
associated with V defined in (2.9) [see also Roe and Pyke (1984), Roe (1986), Harten 
and Osher (1987)]. In the rest of this section, we only deal with W p ; as a result, for 
simplicity of notation, we drop the subscript p. 

Let j be a fixed index. Let V = Vj in (2.8); that is, we linearize the equations 
for the primitive variables at Vj. Since we use a five-point stencil, for —2 < l < 2, 
employing (2.9b), set 

/'-k>/( 2c ) )]“)+! + [!/(2cJ)] p J+ , \ 

W, = LV,„ = pn, - [l/cj] w+ , ; (4.9) 

\ \.P:l (2 c j )] u j+i + [l/(2c|)] Pj+i J 

again, L is the matrix of left eigenvectors at the state Vj. Let tc; be any component 
of W;, say, the first component, for —2 < l < 2. We obtain the final slope for this 
component as follows. Let s+ and s_ be slopes of the linear interpolations: 


s+ = (wi — u>o)/h ; s_ = (u>o — w-i)/h. (4.10a, 6) 

Let p_, po, and p + be the slopes at / = 0 of the three quadratics defined respectively 
by {ia_ 2 ,ta-i,u?o}, {ia-i,it>o,Ufi}, and {tx;o, w 2 }: 

p_ — (w - 2 — 4u;_i + 3u>o)/(2 h), (4.11a) 

po = (ta x - u_x)/(2A), (4.116) 

p + = (-w 2 + 4^1 - 3w 0 )/{2 h). (4.11c) 

Note that the three slopes above are accurate to 0(h 2 ). Next, between / = 0 and 1 = 1, 
the slope s+ does not oscillate, but it approximates the exact slope at l = 0 only to 


0(h). To obtain a slope that is accurate to 0(h 2 ) and does not oscillate, we simply 
use the median function: 

q + = median (s + ,p + ,po), (4.12a) 
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and by symmetry, 


q _ = median (s_,p_,po)- (4.126) 

For further details, see §3 of Huynh (1993). Between the two slopes above, we take, 
roughly speaking, the less steep one, 

q* = xm{q+,q-). (4.12c) 

The slope q m is clearly accurate to 0(h 2 ) because by (4.8a) and (4.2), it lies between 
q_ and q+. It behaves very much like the UNO slope of Harten and Osher (1987): the 
resulting scheme is uniformly second-order accurate and does not create oscillations. 
Because it is (loosely speaking) closest to zero among all slopes that are accurate to 
0(h 2 ), the scheme shares the property of being diffusive with the first-order upwind 
method (slope equals zero) but to a much lesser extent. Note that if the extended 
minmod function in (4.12c) is replaced by the minmod function, then the result is the 
UNO slope. The extended minmod function, however, has the advantage of preserving 
the slope of the smooth side more often near a discontinuity. More importantly, while 
the UNO slope is the final one, the slope (4.12c) is employed only as one of the two 
bounds (the lower limit) for the monotonicity constraint described below. 

Next, we define the other bound (upper limit). The following monotonicity con- 
straint was presented by Van Leer (1974, 1977): for the cell l = 0, if the linear re- 
construction takes values in I [u;_i,u?o, u?i], then it does not create oscillations in the 
solution. This constraint is very simple; however, it causes a loss of accuracy near 
extrema. We present a constraint that preserves uniform second-order accuracy below. 
The argument to the right of / = 0 is carried out first. Van Leer’s constraint can be 
restated as follows. The reconstruction to the right of Z = 0 (xj < x < Xj+ 1 / 2 ) takes 
values in 7[w 0 ,^i]; that is, the final slope q f lies between 0 and 2s+ (see Fig. 4.1). 
It is the limit 0 that causes accuracy to degenerate near extrema. To obtain uniform 
second-order accuracy, we require that the final slope lies between q » and 2s+, where 
q« is defined by (4.12c). The requirement to the left of / = 0 is that q f lies between q, 
and 2 s_. The two requirements together result in the following constraint: the final 
slope lies in the intersection of the two intervals I [<?„, 2s_] and I [<?«, 2s+]. Clearly, one 
end of this intersection interval is q.; the other is 

q * = median (g,,2s_,2s+). (4-13) 
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Figure 4.1. Van Leer’s constraint. 


Thus the constraint takes the form 

9/ € I [?., q*]. 

Having obtained the two Emits, we now define an accurate slope. Using the quartic 
(five-point) formula, set 

q$ = (tu_ 2 — 8u;_i -|- Swi — w 2 )/(12 h). (4.14) 

The above slope is highly accurate; however, near a discontinuity, it may have the 
wrong sign (see Fig. 4.2). We avoid this problem by requiring that q$ Hes between po 
and p m where 

p m = median (p_,p+,p 0 ). (4.15a) 

To bring q$ into the interval I [p 0 , p m ], we once again use the median function: 

q 6 = median (qs,p m ,Po)- (4.156) 

Observe that at the smooth part of the data, q 6 is generally identical to q 5 ; that 
is, p m and po provide plenty of ‘room’ for an accurate slope. Indeed, without loss of 
generality, we may assume Xj = 0, and the Taylor series expansion of the exact solution 
w e takes the form 

w e { x) — sx + dx 2 + ex 3 + 0(h 4 ) 
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x j - 2 x j - 1 x ; x j+l X j + 2 

Figure 4.2. The quartic slope may have the wrong sign. 

[see §3.5 of Huynh (1993)]. Then, by omitting higher-order terms, 

p_ = s - 2 eh 2 , p 0 = s + eh 2 , p+ = s - 2 eh 2 . 

The interval I [p 0 ,p m ] is therefore equal to I [s-2eh 2 ,s + eh 2 ], which contains the exact 
slope s, and the above observation follows. (Also note that instead of the requirement 
that q$ lies between po and p m , we can insist that qs lies between q- and q+, and thus 
save one median calculation. This latter constraint, however, does not provide as much 
‘room’ as the former, i.e., it alters qs more often.) 

Finally, using the median function, we limit q&, 

q f = median (q 6 , q , , q*). (4-16) 

Since q , and q 6 are accurate to 0(h 2 ), the above q f is also accurate to the same order 
by observation (4.4). 

A remark concerning the effect of constraint (4.16) is in order. At the smooth part 
of the data where the slope is nonzero, expression (4.16) yields q / = q& because q. is 
closer to 0 than q 6 , and q * is further from 0. Near an extremum, the interval I [q„ q*] 
may reduce to the point {?*}, and in this case, qf — q * (see Fig. 4.3(a)). It is here 
that our monotonicity constraint preserves second-order accuracy while that of Van 
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x j - 2 x j - 1 x j x j + 1 • r 7+2 x j-2 x j-l X ; + 1 x y+2 

(a) near an extremum (b) near a discontinuity 
Figure 4.3. Effects of the new constraint. 

Leer may not. Near a discontinuity, or where the data change rapidly, the slope qe is 
generally steeper than q *, and the final slope is identical to q *, which is either 2s+ or 
2 s_ (see Fig. 4.3(b)). 

Let the final slopes by the algorithm (4.10-16) for the three components of W; be 
<7^, q^\ and q^p . Let 

Q / = U}\<lP,<lP ) T • (4-17) 

The slopes of the primitive variables are obtained by 

S j = RQ, = RjQ f (4.18) 

where R ; is the matrix formed by the right eigenvectors given in (2.6a). This expression 
completes the slope estimates. 

Note that (4.15) and (4.16) can be replaced by 

?/ = ?*+ minmod (^ 6 - ?», 2 s + — 9,, 2s - - q «). (4.19) 

We also remark that in the case of linear convection, with an appropriate definition 
of monotonicity of the data in an interval (Huynh 1993), it can be shown that the 
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slope defined by (4.10-16) yields a scheme that preserves monotonicity and is uniformly 
second-order accurate (to appear). 

4.2. Smooth regions. The above constraints are constructed so that they have 
little or no effect at the smooth part of the data. For such a region, the slope reduces to 
the quartic formula for the primitive variables. More precisely, suppose the constraints 
(4.15b) and (4.16) have no effect at the jth cell. Then, since q f = q 5 , and multiplication 
by a constant commutes with matrix multiplication, it follows from (4.13) that 

Q/ = [LV;_ 2 - L(8Vj_i) + L(8V i+ i) - LV i+2 ]/(12 h). (4.20) 

Because R and L are the inverse of each other, (4.18) and the above yield 

S; = (Vj-2 - 8V;-! + 8V; +1 - V j+2 )/(l2h), (4.21) 


which is the quartic formula for the primitive variables. 

Next, we present a simple criterion which detects the ‘smooth’ regions where the 
constraints have no effect, and as shown above, the slope is then given by (4.21). For 
a wide class of problems — including all the test problems below — a discontinuity in 
the flow field leads to a discontinuity in the density field. As a result, to separate 
discontinuities from smooth regions, it suffices to test the density field. (For highly 
complex problems, we may need to test more than one field.) Let A jp be the second 
difference of p: 

A|p = Pj- i — 2 pj + pj+ 1. (4.22) 

For each index j, let kj be a flag that is set as follows: if 


1 < < i 

5“ A )p -4’ 


and 


4 tf+iP < 5 

5 - A jp ~ 4’ 


(4.23a, 6) 


then kj = 1, and the flow is considered to be ‘smooth’ in the /th cell; otherwise, kj = 0. 
(We remark that in the case of linear convection, the factors 4/5 and 5/4 are replaced 
by 2/3 and 3/2, and it can be shown that a constraint slightly more general than 
(4.16) has no effect on the original slope. This general constraint was employed for 
cubic interpolation in §3.4 in (Huynh 1993). The factors 4/5 and 5/4 are chosen by 
numerical experiments to account for nonlinearity as well as the fact that only one out 
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of three fields are tested. In fact, the more relaxed factors 2/3 and 3/2 are sufficient 
for all of our test cases except the interacting- blast- wave problem of Woodward and 
Colella (1984).) 

For the purpose of coding, observe that the following statement is equivalent to 
(4.23a): A lies between (4/5 ) Ajp and (5/4) A jp. That is, 

Thus, if we set 

= (!■ a )p - a j- i p) (! a jP - a Up) ’ ( 4 - 24a ) 

and 

h = (f Af, - AJ +1 „) (|A> - A ,j +l p) , (4.246) 

then conditions (4.23) are equivalent to max(ay, bj) < 0. In practice, we allow for some 
noise in the data because they are quickly damped out by the scheme. Conditions 
(4.23) are coded as 

max(a_j, bj) < epj. (4.24c) 

where e is a small number, typically, e = 10~° or 10 -4 . In this paper, all numerical 
tests employ e = 5 x 10 -5 . The corresponding expression for kj is 

k (1 if (4.24c) holds, ( 4 . 25 ) 

1 0 otherwise. 

Roughly speaking, if kj = 1, then the data are ‘smooth’ at Xj, and the slope S ; - is given 
by (4.21). Otherwise, the constraint may take effect, and the estimate for S j is given 
by (4.9-18). 

Notice that at the ‘smooth’ regions (kj — 1), the quaxtic formula (4.21) can be 
replaced by the parabolic one (similar to (4.11b)), which has a three-point stencil; the 
results are generally slightly more diffusive. 

4.3. A slope-steepening technique. The above slope estimate results in a 
scheme that is effective, uniformly second-order accurate, and does not create oscilla- 
tions, but it still slightly smears a contact discontinuity, although to a lesser extent 
than the MUSCL or UNO schemes. The reason for this smearing is that near a dis- 
continuity, the constraints help prevent oscillations by bringing the slope closer to zero 
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when necessary, but they do not steepen the slope to the maximum allowable value, 
which keeps the discontinuity from smearing. In the case of a shock, the characteris- 
tics converge; the shock remains crisp, and no modification is needed. For a contact 
discontinuity, however, the characteristics are parallel; there is no self-steepening mech- 
anism, and smearing takes place. To keep contact discontinuities sharp, Colella and 
Woodward (1984) presented a criterion to detect them, and in that case, a steeper 
reconstruction is employed for their piecewise parabolic method. For other techniques, 
see Harten (1989), Yang (1990), and Mao (1991). Our approach is somewhat different: 
the slope steepening should take effect near a discontinuity but have no effect at the 
smooth part of the data. 

The following simple algorithm can generally resolve a contact discontinuity in four 
cells or less; it is applied only to the linear characteristic, i.e., the second component 
of Wj. For this component, with q~ and q+ given by (4.12), q$ by (4.15b), k a fixed 
constant which typically equals 5, we reset q G to 

q 6 <- sgn (q & ) max(/c|qq - <?_ |, \q&\). (4.26) 

Then we apply (4.16) to obtain the final slope. 

To see the effect of the above expression, observe that at the smooth part of the 
data, since both q- and q+ approximate the exact slope at l = 0 to 0(h 2 ), the quantity 
|< 7 + - q_ | is a small number of order 0(h 2 ). As a result, except for a neighborhood of 
length 0(h 2 ) around an extremum where the original q G is too close to 0, expression 

(4.26) does not alter q & . (We need not pay attention to this neighborhood because the 
constraint (4.16) will make the final slope q f = q..) Loosely speaking, (4.26) has no 
effect at the smooth part of the data. Near a discontinuity, on the other hand, <?_ and 
q + are far apart, and x.\q + — q-\ becomes a large number; consequently, the slope q G by 

(4.26) is considerably steeper than the original q G . 

For readers who are familiar with limiter functions, the following remark explains 

(4.26) from this point of view (see Fig. 4.4). Strictly speaking, the explanation is 
incorrect, but it conveys the idea. Consider the case of q- and q+ both being positive. 
Set k = 5, r = q+/q~, then 


5|g+ ~q-\ = (5|r- 1|)?_. 
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Figure 4.4. Slope steepening explained via limiter function. 

The lines g(r ) = ±5(r — 1) intersect g(r) = (r + l)/2 at r = 9/11 and r = 11/9. The 
lines g(r) = 5(r — 1) and g(r) = 2 intersect at r = 7/5; g(r) = — 5(r — 1) and g(r) = 2r, 
at r = 5/7. Roughly, for r between 9/11 and 11/9, expression (4.26) does not alter 
the original q$, which is represented by g(r ) = (r + l)/2. But for 5/7 < r < 9/11 
or 11/9 < r < 7/5, expression (4.26) steepens q$, and the result can be even more 
compressive than the ‘Superbee’ limiter function (Roe 1985). Also note that since 
both q- and q + are accurate to 0(h 2 ), at most of the smooth regions, r lies between 
9/11 and 11/9, and (4.26) has no effect. Near a discontinuity, it takes effect and 
steepens the slope considerably. The final constraint (4.16) then limits the slope to 
assure that no oscillations are created (r < 5/7 or r > 7/5). 

Certain contact discontinuities, such as those resulting from a shock interacting with 
another shock, are extremely difficult to resolve. Even highly compressive methods like 
Roe’s ‘Superbee’ scheme may smear them. The above technique can be used to capture 
such discontinuities by setting k = 10, or even 20. For such large /c, the slope steepening 
(4.26) may take effect at the smooth part of the data; consequently, it should be applied 
sparingly, for example, only when the strength of the linear characteristic dominates 
that of the other two. This condition is quantified below: for the ith component of 
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w,, 1 < i < 3, set 


<7^ = |ll?2^ — U>o*l + 1^0 * — W -2|- 


(4.27a) 


Then the linear characteristic is said to dominate the other two if 

<tM>2{<tU + <tW). (4.27 6) 

(The factor 2 is found by numerical experiments; the results generally vary only slightly 
if a different factor is chosen.) In this case we also reset kj: 

kj = -2. (4.28) 

As noted at the end of §2.3, all cr (,) have units of density; consequently, a comparison 
among them makes sense. 

We summarize the slope estimate below. 

Algorithm for the slope. Given the conservative variables {U_,}. calculate the prim- 
itive variables { Vj } . and {A jp} by (4.22). For each index j, compute aj and bj by 
(4.24a, b). Next consider the following two cases. 

(a) If (4.24c) is satisfied, the cell is in the ‘smooth’ region; set kj = 1; define Sj by 
the quartic formula (4.21); then move on to the next index. 

(b) Otherwise, set kj = 0, and for — 2 < / < 2, calculate W / by (4.9). For each 
component i, 1 < i < 3, obtain <7^ by (4.27a). If (4.276) holds, i.e., the linear 
characteristic is dominant, reset kj = —2. Next, for each i, 1 < i < 3, calculate 

by (4.10-16); while estimating qf \ however, alter obtaining q& via (4.156), 
if kj = —2, reset qe by (4.26) before carrying out (4.16). Having obtained Q/, 
calculate Sj by (4.18). Then move on to the next index. 

Note that instead of (or in addition to) (4.27b), one can employ a more sophisticated 
criterion to detect a contact discontinuity without much penalty in computing time (at 
least, for scalar computers) because generally, the flow field is mostly ‘smooth’, i.e., 
(4.24c) holds. 

For comparison, we next present the slopes of three popular TVD (total variation 
diminishing) schemes. With s_ and s + defined by (4.10), instead of (4.11-16) the 
minmod scheme (Harten, 1983) is given by 

qf = minmod (s_,s + ); (4-29) 
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Van Leer’s average limiter (1977), 

qj = minmod[(s_ + s + )/2, 2s_,2s+]; (4.30) 

and the ‘Super bee’ scheme (Roe and Baines 1984, Roe 1985), 

qj = minmod[sgn (s_)max(|s_|, |s+|),2s_,2s + ], (4.31) 

We remark that if s_ and s + are of opposite sign, the three expressions above return 0 
for the slope. Expression (4.30) can also be written as 

q 5 = |[sgn (s_) + sgn (s+)] min(|s_ + s+|/2,2|s_|,2|s + |); 

similarly, (4.31) is equivalent to 

qj = |[sgn(s_) + sgn (s+)] min[max(|s_|, |s + |),2|s_|,2|s+|]. 

Compared with the TVD schemes, which are only first-order accurate near extrema, 
the slope estimate for our scheme has the disadvantage of being more involved in terms 
of coding, but it has three distinct advantages: it is both accurate and economical at 
smooth regions, and it gives high resolution at contact discontinuities. 

The author’s SONIC-A scheme, which is a uniformly second-order extension of Van 
Leer’s average limiter, is defined by 

q f = median [(?_ + ? + )/2, ?„?*]. 

The SONIC-B scheme, which extends ‘Superbee’, is given by 

q f = median [sgn (9.) max(|?_|, \q+\), q„ q’}. 

Due to the use of q- and q+ instead of s_ and s + , the above slopes yield uniformly 
second-order accurate schemes. Note that at the ‘smooth’ part of the data, these slopes 
do not simplify easily; in contrast, the slope presented in the algorithm summarized 
above does. 

5. Upwind flux. With the slope {S ; } defined in the previous section, one can 
calculate the time partial derivative by (3.5) and obtain the interface values at time 
£ n+1 / 2 by (3.6a,b). Let j + 1/2 be a fixed interface, and let 

V L = R i (x i+1/2 ,f n+1/2 ), V* = R j+l (x j+l/2 ,t n+1/2 ), (5.1) 
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where the reconstruction values are given in (3.6). From the above left and right states, 
we derive a flux denoted by Ft/ via one of the upwind procedures below. 

In this section, time t is employed only for convenience of derivation; the final flux 
Ft/ is independent of time. Let the conservative variables corresponding to Vt and 
Vr be denoted by Ux and Ur. Consider the conservation laws (2.1) with the initial 
condition at t = 0: 

U = (U L if x < 0, (5.2) 

l Ur otherwise. 

The above problem is called the Riemann problem. The exact solution, which only 
depends on Ux, Ur, and x/t, is known (e.g., Godunov 1959, Van Leer 1979, Holt 
1984, Gottlieb and Groth 1988, Hirsch 1990). It is a constant along each ray (half line) 
x = \ t where the speed A is fixed and t > 0. An approximation Fy for the flux at 
x = 0 (or A = 0) and t > 0, called an upwind flux (or approximate Riemann solver), is 
sufficient for the purpose of calculation. In fact, these approximations are sometimes 
called alternative Riemann solvers because they can be more diffusive than the exact 
solution, and may occasionally produce smoother (more accurate) solutions for (3.3) 
as in the case of a slowly moving shock wave (Roe 1989). Notice that the coordinates 
have been shifted so that (x j+l / 2 , t n+1/2 ) in (5.1) now corresponds to (x = 0,i = 0). 

At an interface in the smooth region, V x and Vr are generally different, but close 
to each other. As a matter of fact, for such an interface, Vr — Vl is an 0(h 2 ) quantity; 
consequently, for simplicity and economy, we may derive F u via the non-conservation 
form (2.3). This method is called the primitive-variable splitting, which can be consid- 
ered as a simplification of Toro’s approximate Riemann solver (1991); it is carried out 
in §5.1. On the other hand, although V L and Vr are close, they are generally different; 
in that case, there is a discontinuity in the initial condition. Theoretically, therefore, 
only the conservation form (2.1) is valid. The simplest model that employs (2.1) is 
perhaps the one-intermediate-state model by Harten, Lax, and Van Leer (1983), which 
resolves the two acoustic waves, but ignores the contact discontinuity. This model 
turns out to be more economical than the primitive- variable splitting; in addition, at 
smooth regions, it performs as well as Roe’s flux-difference splitting with roughly a 35 
percent savings in (total) computing time. 

Near a discontinuity, Vr — Vx is 0(1), and the flux also has an error 0(1). Con- 
sequently, we have to take extra care: if the model is inadequate, this error may 
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propagate throughout the flow field and cause accuracy to degenerate. For such an in- 
terface, we employ a flux-vector splitting derived by linearization and diagonalization. 
It will be shown that this flux-vector splitting, which is conceptually simple, yields a 
result identical to Roe’s flux-difference splitting, which is computationally simple. Our 
derivation also leads to an admissibility correction with no conditional statement and 
tin economical approximation to Osher’s approximate Riemann solver. 

For convenience of notation, with the initial condition (5.2), denote 

AU = Ufl - U L ; (5.3) 

similarly, AF = Fr — Fl, and A pu = Prur — plUl, etc. Notice that the operation A 
is linear, i.e., for any functions /, g, and constants a, /?, 

A (a/ + fig) = aA / + (3Ag. 

5.1. Primitive-variable splitting. Among the subsequent upwind fluxes, this is 
the only one that employs the nonconservation form (2.3). It is perhaps the simplest 
model where the ‘wind’ directions play their roles. As mentioned above, it should 
only be employed at ‘smooth’ interfaces. More precisely, with kj defined by (4.25) and 
(4.28), set 

kj+i/2 = kj + k j+l . (5.4) 

Then this flux is used only when kj +1 / 2 = 2. For a similar method, which employs 
(2.3) with a more sophisticated fixed state for the linearization and a slightly more 
complex calculation involving the left and right states of the contact discontinuity, see 
Toro (1991). 

With the initial condition Yl and Yr, we linearize the Euler equations (2.1) at the 
average state 

V=!(V l +Vk). (5.5) 

Next, equations (2.8a) for the primitive variables are diagonalized as shown in (2.8b). 
The characteristic variables corresponding to the initial condition V l and V r are given 
by (2.9a), 

W p , L = L p V t , and W P , R = L p Yr. (5.6a, b) 
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By (2.8b), each component of these variables is advected with speed V*). Consequently, 
the upwind characteristic variable W Pt u is determined by the ‘wind’ direction or the 
sign of Abb 

w P ,u = 

where, for convenience, the superscript (i) is understood. In both cases, the above 
expression is equivalent to 

w vV = | [l + sgn (A)] w p , L + \ (l - sgn (A)] w p , R 
= \ [w p ,l + w p ,r] - I sgn (A) [w PiR - w PtL ] . 

Denote by sgn (A) the diagonal matrix whose diagonal entries axe sgn(Abl), i — 1,2, 
3; then 

W p ,c = |(W PiL + W P>R ) - | sgn (A)(W p>fl - W P , L ). 

To obtain the corresponding primitive variable V(/, multiply the above equation on 
the left by Rp, and by (5.6a,b), 

Vy — |(Vn + Vr) — |RpSgn(A)L p (AV), (b-7) 

where the operation A is defined in (5.3). Finally, the upwind flux is the flux of the 
above state. 

The following remark is in order. If A = 0, then should we select w p jj — w p> l or 
w P ,u = w p,R in (5.6c)? Since the fluxes resulting from these two quantities are generally 
not identical, the final upwind flux does not depend continuously on the data. At a 
smooth interface, this discrepency is not critical, and the above upwind flux performs 
well; near a discontinuity, however, it may create oscillations in the solutions. By using 
the conservation form (2.1), one can derive upwind fluxes that depend continuously on 
the data; in fact, all of the subsequent methods share this property. 

5.2. One-intermediate-state model. This model, introduced by Harten, Lax, 
and Van Leer (1983), is independent from all other upwind models presented here. 
Instead of using diagonalization, it approximates the solution of the Riemann problem 
(5.2) by three constant states: Ul, U/j, and an intermediate state denoted by U/. Our 


w P ,L if A > 0, 
w PtR otherwise 


(5.6c) 
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presentation below is slightly simpler than the original one. Similar to the primitive- 
variable splitting, this flux is employed at ‘smooth’ interfaces. At ‘nonsmooth’ in- 
terfaces (kj + 1/2 < 2), together with a modification to resolve contact discontinuities 
presented below, it can also be used. 

Our choice of wave speeds is most economical: with the average state V defined in 
(5.5) and c the corresponding speed of sound, let the minimum (subscript /) and the 
maximum (subscript r) wave speeds be 

ai — u — c, a T — u + c. (5.8a, b ) 

The corresponding Mach numbers are 

M = 1, M r = ^=M + l. (5.8 c,d,e) 

c c c 

For a more sophisticated choice of wave speeds, see Einfeldt (1988). Next, the top 
half of the (x,i)-plane is divided into three regions separated by the two half lines 
x = ait and x = a r t (see Fig. 5.1). For the left region, the solution is assumed to 
be the constant state U/,; for the right, U/j; and for the middle region, U/, which is 
determined below. Let r > 0 be a fixed (provisionary only) time step. Consider the 
control volume whose four comers are (a/r, 0), (a;r, r ), (a r r, r), (a r r, 0). By balancing 
the fluxes for (integrating (2.1a) over) this control volume, 

r(a r — a,)U/ = r(a r UR — a/U/,) — r(F/t — Fl). 


As a result, 

TT OrU/i — a/Uf, Fr — Fl . 

U/ = . (5.9) 

a r — ai a r — a/ 

Note that r no longer appears in the above expression. 

We are now ready to derive the upwind flux. For the moment, suppose a; < 0 and 
a r > 0. By balancing the fluxes for the control volume with corners (a/r, 0), (a/r, r), 
(0,r), and (0,0), 

— ra/U/ = -rat U L - t(Fu — Fl). 

As a result, 

Ft 7 = F L + a f (U/-U L ). (5.10a) 
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Figure 5.1. The one-intermediate-state model. 


Similarly, for the control volume with comers (0,0), (0, r), (a r r, r), and (a r r, 0), 


F u = F R + a r (U I -U R ). (5.106) 

Note that the above two F*/ are identical by (5.9). Multiply (5.10a) by a r , (5.10b) by 
— a*, and then add the results, 


Fu = ^-F R + ~^—F L + -^-(Ufl - U L ). (5.10c) 

CLf d\ ai Off Q>1 

Thus the flux F u is a weighted average of the left and right fluxes plus an ‘artificial 
viscosity’ term. Note that the above expression involves only scalar operations and 
no matrix multiplications. Next, if a; > 0 then F u = Fr; similarly if a r < 0 then 
Fy = F r . That is, 


Fu = 


f Fi 

) -=*-] 
\ a r -a[ 


- u t ) 


if 0 < a/, 
if a/ < 0 < a r , 


(5.11a) 


{ F R if a r < 0. 

Observe that if a { = 0, then (5.10c) yields Fy = Fx,; similarly, if a T = 0, it implies 
Fc; = Fr. As a result, the above expression for Ftr depends continuously on the data. 
Using the Mach numbers defined in (5.8c,d,e), 

( F l if 0 < M h 

Fu = { \{-Mi Fr + M r F L ) + lM,M r c(U R - U L ) if Mi < 0 < M r , (5.116) 
1 Fr if Mr < 0. 
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The above expression can be combined into a single formula: 

F u = | [(M r - - Mr) F* + (M r + - M, + )F l ] + \MCM+c [U* - U L ] (5.11c) 
where for any real number z, the negative and positive parts of z are 

z - = min(z,0) = |(z — |z|), z + = max(z,0) = |(z + |z|). (5.12) 

We summarize the algorithm below. 

Algorithm for an upwind flux via the one-intermediate-state model. At a 

‘smooth’ interface where kj+ 1 / 2 = 2, with V i and Vr given by (5.1), first calculate the 
average state V by (5.5); next obtain the Mach numbers M, Mi, and M r by (5.8 c, d, e); 
finally, F u is given by (5.116) or (5.11c). 

Notice that if the maximum and minimum wave speeds are defined as 

ui = — (|u| + c), a r = |u| + c, (5.13a) 

then (5.11a) reduces to 

Ft/ = |(F L + F R ) - |(|u| + c)(U« - U L ), (5.136) 

which is a ‘central difference’ plus an ‘artificial viscosity’ term. This simple formula 
is due to Rusanov (1961) [see also Anderson, Thomas, and Van Leer (1985), Yee 
(1987), and Davis (1988)]. Generally, it saves about 6 percent of total computing 
time compared with (5.11), and it performs well for all of our tests except for the Lax 
problem in §7. 

The one-intermediate-state model (5.11) is economical and produces good results 
at smooth regions. It also handles shocks well; however, it cannot resolve a contact 
discontinuity accurately since the model does not contain a contact. 

We remedy this problem by the following simple technique. At an interface j + 1/2, 
the linear characteristic is said to be dominant if 

kj + 1/2 < 0 (5.14a) 

In this case, we reset pr to pl if u > 0, and reset pl 

pL\ if M < 0, PL*- PR- (5.146, c) 


where kj + 1 / 2 is defined by (5.4). 
to pR if u < 0: 

if M > 0, pr <— 
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In other words, near a contact discontinuity, density is determined by the flow direction. 
Note that velocities and pressures are not altered. The flux F(/ is then calculated by 

(5.5) , (5.8), and (5.11b) or (5.11c). 

With the above modification, the model works well for all the standard tests except 
the Mach 3 expansion shock in §7; i.e., it may admit non-physical solutions (also known 
as solutions that violate the entropy condition). Further modifications can be made 
so that non-physical solutions are excluded; e.g., at a ‘nonsmooth’ index j + 1/2, i.e., 
kj+ 1/2 < 1> instead of (5.5), a/ and a T can be defined by 

ai = min(u£ — ci,up — cp), a r = max(ui, + c^up + cp). 

Such calculations involve two speeds of sound and take more computing time than 

(5.5) . We will not elaborate on different definitions for a/ and a r any further; instead, 
we present below an estimate for F jj that has two intermediate states, and with an 
admissibility condition, it excludes non-physical solutions. 

5.3. Flux-vector splitting. The flux-vector splitting developed by Sanders and 
Prendergast (1974) and Steger and Warming (1981) is conceptually simple; however, 
the resulting scheme smears a contact discontinuity and may become inaccurate at 
sonic points [see also Hirsch (1990)]. On the other hand, Roe’s flux-difference splitting 
(1981, 1986) is conceptually more involved, but it is highly accurate. The models for 
these two splittings are very different: for the former, the interaction between the left 
and right states is accomplished through mixing of pseudo-particles (the Boltzman 
approach); for the latter, through finite-amplitude waves (the Riemann approach) [see 
Haxten, Lax, and Van Leer (1983)]. 

We present below a flux- vector splitting that turns out to be identical to Roe’s flux- 
difference splitting. The advantage of this presentation is its conceptual simplicity. The 
method, in fact, is very similar to that of §5.1, but instead of splitting the primitive 
variables, we split the fluxes. The two key techniques are, again, linearization and 
diagonalization. What is crucial, however, is that the fixed state for the linearization 
must be derived carefully in order that the solution depends continuously on the data. 
Our derivation of this state, which is Roe’s ‘square root averaging’ state, follows the 
remark concerning M. Brio in Roe and Pyke (1984). 
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As expressed in (2.10), after a Linearization, the Euler equation can be written as 


<9U - d\J 

dt c dx 


— 0 ; 


here A c = (<9F/5U)(U), and U is yet to be determined. With the initial condition 
( 5 . 2 ), let U be a state that satisfies 


AF = A C AU 


(5.15) 


where A is defined in (5.3). (As will be shown later, the above condition implies the 
upwind flux depends continuously on the data.) Loosely speaking, we wish to find a 
solution for the mean value theorem (which says a solution exists in the scalar case). 
Using expression (2.22) for A c , 

f A pu \ / 0 1 0 \ / A p \ 

A {pu 2 + p) = (7 - 3)tt 2 /2 _ (3 - 7 )u _ 7-1 A pu . (5.16) 

^ A(e + p)u J \ (7 — l)u 3 /2 — uH —(7 — l)u 2 + H 7 u ) \ Ae j 

The first component of the above vector equation yields A pu = A pu, which is useless; 
as a result, the solution of the above system cannot be unique. Next, using expression 
( 2 . 2 b) for e, the second component can be simplified to 

u 2 Ap — 2uA(pu) + A (pu 2 ) = 0. (5.17) 


The two solutions of the above quadratic in u axe 

A (pu) ± [{A(pu )} 2 - ApA(pu 2 )] 1/2 


u — 


Ap 


For simplicity of notation, set 


r L = \fph-, r R = yfpR, 


(5.18) 


then, by expanding A via (5.3), 


n/OL ± rRUR 
r L i r R 


(5.19) 


The solution corresponding to the negative sign is meaningless when pl = pn , and is 
therefore discarded. The one with the positive sign always makes sense and yields 


tlul + r R u R 

u = 

r L + r R 


(5.20) 
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The above expression is 

Roe’s well-known ‘square root averaging’. 

For the purpose of 

coding and for later use, denote 



n rL n rR 

h>L — . 1 KK — . 5 

tl + tr rL -t rR 

(5.21a, b) 

and 

P* — V PL PR- 

(5.22) 

Then 

0L = - — > ^r — I-^l-, 

pL + P* 

(5.23a, b ) 

and 

u = Plul + Prur', 

(5.24) 


that is, u is a weighted average of ul and ur. 

Consider now the third component of (5.16). Since u is already determined above, 
this component yields the total enthalpy H. Indeed, by definition (2.17) for H, 


e + p — pH, (5.25a) 


and by (2.2b), 


e = -pH + — — — pti 2 . 
7 2 7 


(5.255) 


Using the above expressions, (5.17), and the linearity of A, the third component can 
be simplified to 

A (puH) = —uHAp + HA(pu) + uA(pH). 


As a consequence, 


H = 


A(puH) — uA(pH) 


A(pu) — uAp 

By expanding A in the above expression and factoring out ul ) ■> 

tlHl + trHr 


H = 


rL + r R 


that is, with and &r by (5.23a,b), 

H = p L H L + PrHr 


(5.26) 


(5.27) 


(5.28) 


(note the similarity between (5.28) and (5.24)). 
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The argument (5.16-28) yields the following important fact: equation (5.15) is 
satisfied with u by (5.24), H as above, and p arbitrary. 

Next, we split the fluxes. As noticed after (2.20), we do not need p for the purpose 
of diagonalization because it does not appear in the expressions of A c , A c , L c , and R c . 
Multiplying (2.1) on the left by L c , and since L c is a constant matrix, 


dL c U dL c F 

+ = 0 . 


dt ' dx ^ 5 ‘ 29a) 

Let the characteristic variable be W = L C U (for convenience, we omit the subscript c 
in W), and the characteristic flux, G = L C F. Then, the above is equivalent to 

dW dG 


dt + dx °‘ 
Corresponding to the initial condition (5.2), 

W L = L c U l , W fi = L c Ufi, 

and more importantly, 


(5.296) 

(5.30) 

(5.31) 


Gl = L c Fl,. Gh = L c Fi*. 

Since L c is a constant matrix, (5.30) and (5.31) respectively imply 

AW = L C AU, AG = L C AF. (5.32a, 6) 

We now make use of the crucial expression (5.15): AF = A C AU. Substitute this into 
(5.32b), 

AG = L C A C (R C L C )AU. 

By (2.15) and (5.32a), 

AG = AAW. (5.32c) 

Let be the ith component of G, then the upwind characteristic flux G^- at x = 0 
for (5.29b) is determined by the sign of 


Su = \ 9L if t . A -°’ 

( gn otherwise 


(5.33) 


where, once again, the superscript (?) is omitted. Expression (5.32c) shows that if 
A = 0, then gi = gn, and it does not matter whether the left or the right flux is 
selected in (5.33). Finally, the upwind flux is 

F u = R C G[/. (5.34) 
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In summary, let p, be defined by (5.22); /?£,, 0r by (5.23); u, (5.24); H, (5.28); L c , 
(2.18); and R c , (2.16). Then the upwind flux is given by (5.31), (5.33), and (5.34). 

Similar to the derivation of (5.7), the above algorithm can be expressed with no 
conditional statement. Indeed, by considering the two cases of (5.33), 


gu 


i + 


sgn (A)]^l + 5 [l-sgn(A)]$fc. 


As a result, switching to vector notation, 


(5.35) 


Gy = + Gr) — |sgn — Gl)- 

Denote by |A| the diagonal matrix whose diagonal entries are |A^|. Then from (5.32c), 
and since [sgn (A)] A = |A|, 


G(/ = \(G l + Gr) - ||A|AW. (5.36) 

By (5.32a), 

Gy = 5(Gl + Gr) — | |A| L C AU. (5.37) 

Multiplying the above expression on the left by R c and employing definition (5.31), 
one obtains 

F v = |(F L + F r) - |R C |A| L C AU. (5.38) 

Thus, in the flux- vector splitting algorithm summarized after (5.34), expressions (5.31), 

(5.33) , and (5.34) can be replaced by (5.38). Expression (5.38) is very similar to Roe’s 
flux-difference splitting formula; in fact, in the next subsection, we will show that they 
are identical. 

Note that in (5.31) and (5.33), we split the left and right fluxes to obtain the upwind 
flux; in (5.38), however, it is AU, the difference of the conservative variables, that is 
split. As shown above, the two splittings are identical provided (5.15) is satisfied. If a 
different state is chosen for the linearization, e.g., the average state U defined in (5.5), 
the splitting (5.33) no longer yields a result identical to (5.38). 

With the tilde state given by (5.15), expression (5.33) yields a flux that depends 
continuously on the data. If a different state is chosen for the linearization, the flux via 

(5.33) may not depend continuously on the data any longer. In contrast, expression 
(5.38) always results in a flux that depends continuously on the data no matter which 
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state we pick for the linearization. For this reason, (5.38) is preferred. Conceptually, 
however, (5.33) is simpler. 

Also note the similarity between the primitive-variable splitting (5.7) and the flux 
splitting (5.38). What is crucial, however, is the differences between the two: to the 
advantage of the flux splitting, the sign function in (5.7), which is discontinuous, is 
replaced by the absolute value function in (5.38), which is continuous; to its disad- 
vantage, the matrices in (5.7) associated with the primitive variables are sparse, while 
those in (5.38) associated with the fluxes are not. 

Observe that the key difference between the flux-vector splitting of (5.33) and that 
of Steger and Warming is that in the latter, the left flux is split by the matrix L Cy i at 
the left state, and similarly, the right flux, by the matrix at the right state; whereas 
in the current method, both the left and right fluxes are split simultaneously by the 
matrix L c at the tilde state. 


5.4. Flux-difFerence splitting. In this subsection, we show that the above 
method yields a result identical to Roe’s flux- difference splitting. To prove this, we 
first calculate p by employing the Jacobian of the transformation between the primitive 
and conservative variables. In addition to (5.15), we require that the tilde state satisfies 


AU = MAY, 


(5.39) 


or, with M expressed by (2.11b), 


Ap\ 

A pu I = 
Ae j 


1 

u 

W 2 


0 

P 

pu 


0 

0 

1/(7 - 1) 



(5.40) 


Notice that due to the use of M, our derivation of the tilde state is conceptually simpler 
than that of Roe and Pyke (1984), where the eigenvectors are employed to decompose 
AF and AU. 

The first component of the above vector equation yields A p = A p, which is again 
useless. The second and third yield, respectively, 


A (pu) = uAp + pAu, 

(5.41) 

A (pu 2 ) = u 2 Ap + 2puAu. 

(5.42) 
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Note that these two equations are identical to (3.9-10) of Roe and Pyke (1984). Cal- 
culating p from (5.41) and substituting it into (5.42), one obtains 

u 2 Ap — 2uA(pu) + A (pu 2 ) = 0. 

Strikingly, this quadratic is identical to (5.17), and the solution for u is given in (5.20). 
Substitute u into (5.41), one obtains 

p = y/PLpR ■ (5.43) 

The argument (5.40-43) yields the following important fact: equation (5.39) is 
satisfied with p as above, u by (5.24), and H arbitrary. 

As a result of the statement above and that after (5.28), the following crucial 
conclusion can be drawn. There exists a state V which satisfies the two equations 

AF = A C AU and AU = MAV (5.44a, b ) 

simultaneously, and this state is unique: it is defined by (5.43), (5.24), and (5.28). 

We make use of (5.44b) to simplify (5.38). Multiplying (5.44b) on the left by L c 
and employing (5.32a) and (2.14b), 

AW = L C AU = LpAV. (5.45) 

Substituting this into (5.38), we obtain 

Ft T = |(F L + Ffl) - §R C | A| L P AV. (5.46) 

The above is identical to expression (28) of Roe (1986), which is the flux-difference 
splitting formula. To clarify the differences in notation, note that the components of 
LpAV are denoted by a k by Roe; and the columns of R c , by e^. Since the flux-vector 
and flux-difference splittings yield identical results, when a distinction between them 
is not necessary, they are called the flux-splitting method. 

Observe that the advantage of the flux-difference splitting (5.46) over the flux split- 
ting (5.38) is that the matrix L p is sparse while L c is not. 

Also note the similarity between the flux-difference splitting and the primitive- 
variable splitting (5.7). For comparison, suppose we replace the tilde state in (5.46) 
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by the average state (5.5). Then, compared to the primitive- variable splitting, the 
flux-difference splitting has the disadvantage of costing slightly more since Rp is sparse 
while R c is not; on the other hand, its advantage is that the resulting flux depends 
continuously on the data, and consequently, it is more robust. 

We remark that with the tilde state defined by (5.24), (5.28), and (5.43), one can 
derive p: 

P = PrPl + PlPr + — ^l^kp(Au) 2 . 

Z7 

From the definition of the tilde state, it is not obvious that p is strictly positive. But 
from the above expression, and since p by (5.43), /?l, and 0r by (5.23) are all strictly 
positive, so is p. 

As discussed by Roe (1981, 1986, 1989), the flux-splitting method may admit non- 
physical solutions. This problem can be corrected by an admissibility condition (Roe 
and Pyke 1984) or by adding numerical diffusion (Harten 1983, Harten and Hyman 
1985). One can also derive an upwind flux that satisfies the entropy condition a 
priori (Engquist and Osher 1981, Osher and Solomon 1982, Osher 1984, Osher and 
Chakravarty 1984). For a recent treatment, see Roe (1992). Two types of admissibil- 
ity conditions that are somewhat different from the one in Roe and Pyke (1984) are 
shown below. Our presentation also establishes a connection among these different 
approaches: they all add numerical diffusion by increasing the magnitude of the wave 
speeds. 

5.5. Admissibility conditions for flux-vector splitting. In this subsection, 
we present a conceptually simple admissibility condition that excludes non-physical 
solutions. Due to its simplicity, it may easily extend to other systems of equations. 

Again, the superscript (i) is understood. Suppose A l and Xr (yet to be determined) 
are respectively the speeds associated with the characteristic quantities {wl^l] and 
defined in (5.30-31). Then, for the flux-vector splitting, we have effectively 
assumed that Xl = Xr = A. In other words, the *th component of equation (5.29b) is 
a convection equation with speed A. Since (5.29b) is nonlinear, if we can approximate 
X L and Xr more accurately, in addition to the case of a convection, we can have either 
that of a shock where Xl > X > Xr, or a fan where Xr < X < Xr. 

In the case of a shock (see Fig. 5.2(a)), expression (5.33) is still appropriate for 
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(a) shock (b) fan 

Figure 5.2. Characteristic fluxes for a shock and a fan. 

the following reasons. The upper half of the (x, t)-plane is divided into two regions by 
the shock line x = At, and the characteristic flux is a constant for each region, taking 
the value gL to the left of the shock, and the value gp to the right. Expression (5.33) 
determines which of these two values is chosen for gu, the characteristic flux at x — 0 
(x = At, A = 0, and t > 0). 

In the case of a fan (see Fig. 5.2(b)), however, (5.33) may no longer be appropriate 
for the following reasons. The upper half of the (x, t)-plane is divided into three regions 
in this case: to the left of x = A x,t, the flux is gL\ to the right of x = Xp f, it is gp, and 
for the region consisting of the rays x = At where A l < A < Xp, the flux g{X) changes 
smoothly from gL to gp. If the fan does not spread past x — 0 (or A = 0), then clearly, 
(5.33) is still valid. But if the fan spreads past x = 0, (5.33) still yields either g i, or gp , 
which is a poor approximation for gu and may result in non-physical solutions. This 
situation is depicted further in Fig. 5.3, where the characteristic flux g as a function of 
the characteristic variable w is plotted. 

We estimate A l and A# by using the tilde state as a pivot, and approximate the 
fan by a (piecewise) linear function below. First, denote the vector of characteristic 
variables and fluxes at the tilde state by W and G: 

W = L C U = L P (M-'U), G = L C F = L P (M _ 1 F). (5.47) 

By carrying out the multiplications in the brackets using (2.11c) and then the multi- 
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Figure 5.3. A fan spreading past A = 0 in {w^g) -space. 


plication by L p using (2.6c), the above expressions yield 


W = — 

7 


1/2 \ 
7 - 1 
1/2 


G = — 

7 


f {u — c )/2 

«(7 ~ 1) 
\{u + c)/2 t 


(5.48a, b ) 


(We remark that G r and Gr defined in (5.31) can be simplified in a similar manner.) 
Note that G can also be derived by observing that the flux F, as a function of the 
conservative variable U, is homogenous of degree one: F = A C U; as a result, G = AW. 
In addition, notice that each component g of the vector G has units of momentum, 
while each component w of W has units of density and is strictly positive. Next, let 
the average speeds be 


Az, 


9L_Zl 

WR — u?' 


Xr — 


9R ~9 

Wr — w 


We can now estimate and A^: 


Xr = 2 A r — A, A r = 2Xr — A. 


(5.49a, b) 


(5.50 a, 6 ) 


If A > 0 and Xr < 0, then we have a fan that spreads past x = 0. By assuming that the 
flux g is a linear function of the speed A between (Xr^r) and (A, < 7 ), we can estimate 
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gi 7, the flux corresponding to the speed 0: 

gu - 9l 


9-9l 


—X l A — Al 

Thus, for this case, we reset gu defined in (5.33) by 

A gL — A l9 


gu 


Similarly, if A < 0 and Xr > 0, then 


gu 


x~x L 


— XgR + A Rg 


(5.51a) 


(5.516) 


—A + A r 

By considering the limiting case of Xi — > 0 and that of A r — ► 0, one can easily verify 
that gu defined by (5.33) and (5.51a, b) depends continuously on the data. 

Notice that we need to carry out the above admissibility correction only for the 
nonlinear characteristics; for the linear one, (5.33) is sufficient (more on this in the 
next subsection). In fact, if u > 0, then the first component g^ requires this cor- 
rection, but not the third since the left and right speeds corresponding to a nonlinear 
characteristic do not spread beyond the contact speed u (approximately anyway; for 
the exact solution, this is always the case). Similarly if u < 0, then the correction is 
carried out only for gjfK 

Also note that in general, X L and A r given by (5.50) may lie on the same side 
relative to A. However, if the left and right states differ only by either the first or the 
third characteristic, then, for all practical situations, Xr and A r lie on different sides 
relative to A; i.e., we have either a shock or a fan. 


5.6. Admissibility conditions for flux-difference splitting. Since the role 
of the admissibility condition is to break up the solution when it is not physical, a 
simple model may serve the purpose. In this subsection, we make use of the simple 
wave model employed by Roe and Pyke (1984) to approximate X L and Xr, and then, 
derive an admissibility condition for the flux-difference splitting that involves no con- 
ditional statement. Moreover, formula (5.67) below shows that this condition produces 
additional ‘numerical viscosity’ when the fan spreads past x = 0, and thus, (5.67) 
provides an explanation for the connection between the approach of admissibility con- 
dition (Roe and Pyke 1984) and that of numerical diffusion (Harten 1983, Harten and 
Hyman 1985). 
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_ A 

Let V = (/>, u, p) T be a fixed state, which plays the role of V in (2.9). Let V + 6V = 
(p+8p, u+8u, p+8p) T be a nearby state. For each ith component, we will estimate 6Ab), 
the change in the characteristic speed, in terms of 8w^\ the change in the characteristic 
variable. By linearizing the Euler equations at V, (2.9) implies 


8W = L P 8Y = 


~[p/ (2c)] 8u + [1/ (2c 2 )] 8p 
8p — [1/c 2 ] 8p 
[pI (2c)] 8u + [1/ (2c 2 )] 8p 


(5.52) 


(See also expression (5.45).) 

Next, we assume that the two states differ from each other only by the first char- 
acteristic; that is, 

8w^ = 8p — [1/c 2 ] 8p = 0, (5.53a) 


and 

<W 3) = [p/(2c)\ 8u + [1/ (2c 2 )] 8p = 0. (5.536) 

All 8 quantities are expressed in terms of 8p below. Equations (5.53a, b) are equivalent 
to 

8p — c 2 8p, 8u = —(c/p)8p. (5.54 a, 6) 

As a consequence of (5.52) and (5.54), 


8w = 8p. 


(5.55) 


Now, we calculate 8c. By definition (2.4) for the speed of sound, 

2c6c = 

p 2 

Using (5.54a) and again (2.4), the above implies 

(7 - 1 )c8p 


8c -- 


2 P 


(5.56) 


The speed corresponding to the first characteristic of V is = u — c; that of V + 8~V 
is A^U + 6A(U = u -| - 8u — (c + 8c). Consequently, by (5.54b) and (5.56), 

6\m = 6u 

2p 
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And by (5.55), 


jjdl = _ (7 + lHt 0 < 1 > (5.57) 

2p 

A similar argument for the third characteristic with the assumption that Sw (1) = 0 
and <W 2 ! = 0 yields 

fA<»> = h ± 1 )^ u,(3> . (5.58) 

2 ? 

For the second characteristic, = 0 and Sw^ = 0 imply that Su = 6p = 0. As 
a result, £ A^ = 6u = 0. For this reason, the second component needs no modification. 

We are now ready for the admissibility correction. With the initial condition V x, 
and Vfl, let the tilde state be defined by (5.43), (5.24), and (5.28). Let AW = L P AV, 
and for the purpose of derivation, let G u be defined by (5.36). 

We carry out the argument for the first and third characteristics below. Again, the 
superscript (i) is omitted unless the corresponding expressions for the two characteris- 
tics are different. Equation (5.36) implies 

gu = \(9L + 9R) - ( 5 - 59 ) 

For the purpose of estimating A L and X R , we assume that w = (w L + w R )j 2; therefore, 

wr — w = w — wl — | Au;. 

Next, we employ either estimate (5.57) or (5.58): with V playing the role of the fixed 
state V; V L and Vr, that of V + 6V; in the case of the first characteristic, set 

AA = AAW = - (7 + 1 >^"' -, (5.60a) 

2p 

or if it is the third one that we are dealing with, set 

AA = AA (3) = (7 + 1 0 } ' A - U --. (5.605) 

2p 

Then 

A l = A-|AA, Afl = A + |AA. (5.61a, b) 

Note that if AA < 0, then X L > X > X R , the characteristic speeds converge, and we 
have a shock. On the other hand, if AA > 0, then Ax, < A < X R , the characteristic 
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speeds diverge, and we have a fan. Furthermore, for the first characteristic, AA and 
Aw have opposite sign by (5.60a); for the third, the same sign by (5.60b). 

Denote the speed 0 by A 0 and the corresponding characteristic variable and flux by 
wo and go, respectively. (Thus go = gu •) Suppose now we are in the situation of a fan 
that spreads past x = 0; i.e., A^ < 0 < Xr. In addition, suppose A > 0. It follows from 
(5.33) that expression (5.59) yields gL- Since the correct choice for gu is go, we should 
add the quantity g 0 — gL to the right hand side of (5.59). 

To keep the formulas simple, we estimate <70 — gL in terms of Aw and the various 
speeds below. The following estimates are independent of the assumption of a fan 
described above. As assumed in (5.60a-61), the speed A is a linear function of w ; 
therefore, 

wo — wl _ A 0 — A l _ —A l 
A w ~ AA ~ ~AA~' 

Moreover, wl corresponds to the speed A^, and wq to 0; consequently, the average 
speed between w 0 and wr can be approximated by A^/2: 

go — gL ^ A l 
w 0 — w L 2 ' 

As a result of the above two expressions, 


A similar argument yields 


go - gL 


go - gn 


A. 

2AA 


Aw. 


X 2 R A 

-Aw. 


2AA 


(5.62a) 


(5.626) 


Now, we turn to the situation of a fan that spreads past x = 0. There are two 
subcases: for A > 0, set 

77 = (5.63a) 

for A < 0, set 


T] = 


AX' 


(5.636) 


The corrected upwind flux is given by adding go — gL to the right hand side of (5.59) 
in the former case, and by adding go — gR in the latter. As a result, by (5.62-63), 


gu = \ G 9l + 9r) - |(|A| + g)Aw. 


(5.64) 
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Next, we derive a formula for rj that works for all cases. Instead of A l and \r, 
consider the quantity 

|A|-|AA. 

In the case of a shock, i.e., AA < 0, the above expression is positive. For a fan that 
does not spread past x = 0, it is also positive. For a fan that spreads past x — 0, it 
is negative; it equals A L if A > 0; but -A r if A < 0. Only when the fan spreads past 
x = 0, does the flux need a correction. Thus, with 

A;v = {|A| — |AA} _ , (5.65) 


where z is the negative part of z defined in (5.12), set 

A % 


T) = 


|AA| + 10 


—15 1 


(5.66) 


then the corrected flux is given by (5.64). 

Let H be the diagonal matrix whose diagonal entries are ?7 (1) , 0, and r/ (3) . By 
multiplying the matrix expression corresponding to (5.64) on the left by R c , one obtains 


F a = \ (F t + F*) - JR, (|A| + H) AW. (5.67) 


Before summarizing the algorithm, note that the quantities 7-I, l/(q-l), 7/ (7—1)1 
and (7 — l)/7 are often used and should be stored. 

Algorithm for flux-difference splitting with admissibility correction. At a 

‘nonsmooth ’ interface, i.e., ^+1/2 < 2 where kj+ 1/2 is defined by (5.4), with V+ and 
\ R given by (5.1), calculate: and eji via (2.2b); Hl and Hr, (2.17); p, (5.43); /?l 

and 0 R , (5.23); u, (5.24); H, (5.28); c 2 , c, (2.17); AV, (5.3). With L p by (2.6 c), obtain 
AW = LpAV. Next, for i = 1 and i = 3, calculate A (t ) via (2.5); AA^, (5.60a); AA (3) , 
(5.606); a{J, (5.65); 17W, (5.66). Finally, with R c by (2.16), F L and F r, (2.2c), the 
upwind flux is given by (5.67). 

Observe that expression (5.66) can be replaced by a simpler one via an additional 
approximation as shown below. Rewriting (5.66) as rj = [ Ajv |(|A jv|/|AA|), then by 
(5.65), the quantity in the parenthesis is bounded between 0 and 1/2. We can replace 
this quantity by its upper bound 1/2: 


7 = ||V|- 


(5.68a) 
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The above expression always yields a more diffusive flux than (5.66). Instead of the 
factor 1/2, we can use 1/4, which is the average of 0 and 1/2, 

? = i|A*|. (5.686) 

For the test problems below, the results by these two expressions axe essentially iden- 
tical to those by (5.66); nevertheless, only (5.66) is employed here. 


5.7. Simple-wave upwind flux. We can now derive an economical approximation 
to Osher’s approximate Riemann solver (Engquist and Osher 1981, Osher and Solomon 
1982). First, in our context, instead of (5.59), this flux is expressed as 

9u = r {gi + 9r) ~ ~ / \M dw. (5.69) 

Z Z Jwl 

Note the key difference between Osher’s method and the above expression: for Osher’s 
flux, the integral follows a certain path in the space of state variables, and determining 
this path is quite involved; here, the path is already determined by the linearization 
and diagonalization. Observe that symbolically, 

da f [ WR 

— = A, /A dw = g, / A dw=g R -gL- 

dw J J iv L 

As a consequence, if Ax, >0 and A r > 0, then since A is assumed to be linear in u?, 
A(u>) ^ 0 for all u) between to t , and to^j^ and expression recovers gL) similar ly^ if 

Al < 0 and < 0, it recovers gn. For a fan that spreads past x = 0, it yields 


l 1 o 1 r w R 

9u — - ( 9i + Qr) + - / A dw — - Xdw = go 

Z Z Jwl Z J wo 


(5.70) 


For these three cases, the result by (5.69) is the same as that by the flux-difference 
splitting with admissibility condition in the previous subsection. 

However, if we have a shock with left and right speeds spreading past x = 0, i.e., 


Xi > 0 and \r < 0, 


(5.71) 


then (5.69) yields 

1 1 f w ° 1 f WR 

9u = t:{9l + 9r) ~ o / A dw + - Xdw = g L + g R - g 0 . (5.72) 

Z Z Jwl Z J XU® 
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Next, in addition to (5.71), we assume that A > 0. Since (5.59) yields g L , we must add 
to its right hand side the quantity gji — go to obtain (5.69). But this quantity is given 
by (5.62b), and since AA < 0, 


9R- 9 o 


2|AA| 


Au;. 


(5.73) 


Thus, with 


rj = 


|AA|' 


(5.74a) 


the upwind flux resulting from (5.69) is given by (5.64). The case A < 0 is similar 
except 


g = 


2L 

|AA|‘ 


(5.746) 


An expression that works for all cases can be derived by assuming that we always 
have a fan in (5.65): A A/2 is replaced by its absolute value, 


A;v = {|A| - ||AA|}-. (5.75) 

Clearly, for a fan , the above expression reduces to (5.65) since AA > 0. For a shock 
that does not spread past x = 0, it yields 0. For a shock spreading past x = 0, there 
axe two subcases: if A > 0, it yields A#; if A < 0, it yields —\l- 

Thus, with X N defined by (5.75) instead of (5.65), the algorithm for the flux- 
difference splitting with admissibility correction in the previous subsection yields an 
approximation to Osher’s approximate Riemann solver. This method is named the 
simple- wave upwind flux. 

Observe that the argument (5.72-74) can also be carried out by switching the limits 
of integration in (5.72) so that the variable A corresponding to a shock is transformed 
into that of a fan. 

Note that the flux-difference splitting turns a fan into a shock, whereas the algo- 
rithm in this subsection does the exact opposite by turning a shock into a fan. Also 
notice that to save computing time, one can replace the tilde state by the average state 
(5.5). 


Flow with area variation. The formulation presented in the previous sections 
extends easily to the case of flows with source terms. As an example, consider the 
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unsteady one-dimensional flow in a duct of variable cross-sectional area A(x), which is 
described by 

(AU), + (AF) r = B, (6.1a) 

where U and F are defined in (2.1), and the source term is 

B = (0,pdA/dx,0) T . (6.16) 

This problem has been studied by several authors; see for example Roe (1986) and the 
references given there. In the same paper, Roe presented a modified Riemann solver 
that takes into account the area variation. Using Hancock’s half-step formulation pre- 
sented in section 3, the area variation is accounted for simply in the partial differential 
equations for the primitive variables derived from (6.1). First, (6.1) implies 

U t + F x + = 0, (6.2a) 

where 

C = F — (0, p, 0) T = (pu, pu 2 , (e + p)u) T . (6.26) 

Next, with V and A P defined by (2.3), the above leads to 

V t + A P V x + = 0, (6.3a) 

where 

K = (pu, 0,7pu) T . (6.36) 

The first component of (6.3a) follows easily from that of (6.2a); a similar statement 
holds for the second component, which employs the first; and the derivation of the 
third component employs both the first and the second. 

For each jth cell, the calculation of Sj, which approximates (V x )j, is the same 
as that of section 4; however, the quantity Tj, which approximates (V t )y and was 
estimated by (3.5), is now obtained by using (6.3): 

Tj = — (A,)j S,- — (^~^j . K,. (6.4) 

This is where the source term plays its role. The rest of the calculations are the same 
as those for the Euler equations (2.1): the solution at time f n+1 / 2 is given by (3.6), 
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and the flux at each interface is obtained by one of the methods in §5 without any 
modification because time is frozen at t” +1/2 , and the source term contributes nothing 
here. Finally, using the midpoint rule in a manner similar to the derivation of (3.3b), 
expression (6.1) implies 


Aj Uf 1 = A, U“ + £(A,-- 1/2 F£$ - A i+1/2 + Z," +1/2 ), 


-,n+l/2 


n+l/2> 


(6.5a) 


where 

Zf 1/2 = (o, p ” +i/ 2 [A j+1/2 - A,-.,/,]. o) T , (6-56) 

f>r /2 = Pi + Mp*);, («- 5c ) 

and (pt)j is given by the third component of (6.4). The above expression completes 
the algorithm for the interior points. Note that in the steady state case, one can set 

n+l/2 

Pj =Pj- 

Our boundary conditions are the standard ones. At the inflow boundary, if the flow 
is supersonic, we specify all three variables p, it, and p; if it is subsonic, we specify two 
and extrapolate the third from the interior. For the outflow boundary, if the flow is 
supersonic, we extrapolate all three variables] if it is subsonic, we extrapolate two and 
specify one, say, pressure. For simplicity, and since the flow fields in our numerical 
examples are essentially constant near the boundary, we employ a piecewise constant 
reconstruction for the boundary cells. For the cell right next to the boundary ones, 
the slope is estimated by the parabolic (three-point) formula. One can also employ 
appropriate one-sided interpolations; the results are essentially the same. 


7. Numerical experiments. Unless otherwise stated, in all of our numerical ex- 
amples, 7 = 1.4; the slope steepening coefficient k is set equal to 5; at a ‘smooth’ inter- 
face where fcj+1/2 = 2, the one-intermediate-state model is used; otherwise (fcj+1/2 < 2), 
the flux-difference splitting with admissibility condition is employed; and finally, the 
CFL number is 0.8. Note that at a ‘smooth’ interface, the computing effort is only 
slightly more than that of the central difference scheme with artificial viscosity (5.13b) 
[see also Jameson, Schmidt, and Turkel (1981), and Harabetian and Pego (1993)]. Gen- 
erally, most of the interfaces are ‘smooth’. Near a discontinuity, i.e., at a ‘nonsmooth’ 
interface, the computing effort is slightly more than that of the UNO scheme. Also 
note that for convenience, the one-intermediate-state model is said to be employed 
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everywhere if at a ‘nonsmooth’ interface (kj+ 1/2 < 2), this model with the modification 
(5.14) is used. 

For shocktube problems, the initial time steps corresponding to a fixed CFL condi- 
tion are often too big because the flow has not developed; e.g., for the Riemann problem 
used by Sod (1978), with CFL number 0.8, the first three time steps are 0.0068, 0.0040, 
and 0.0039, while the final time step is 0.0036. Because the solution is self-similar, all 
time steps are equal to 0.0036 in the exact solution. The initial time steps that are too 
big can occasionally cause some oscillations, especially for CFL number 0.9 and above. 
To avoid this problem, one can make a test run to find the appropriate time step, and 
then fix the time step rather than the CFL condition. This modification, however, is 
not appropriate in the case of Woodward and Colella’s blast wave problem, where the 
time steps vary as time goes by. The following simple modification works well for all 
test problems below: let the kth. time step corresponding to a fixed CFL condition be 
denoted by At; for the first 9 time steps, (1 < k < 9), we reset At by 


At 


At 



(10 - k) 2 ' 
100 


(7.1) 


Generally, to arrive at the same final time, a test case with the above modification 
takes three time steps more than one without. For the example mentioned above, the 
first three time steps are 0.0013, 0.0020, and 0.0023. Expression (7.1) is employed in 
all of the following problems. 

The first test, used by Sod (1978), is the Riemann problem (5.2) with the initial 
data 

(pl,ul,Pl) = (1,0, 1), (pR,u R ,p R ) = (0.125,0,0.1) (7.2) 

where the subscript L represents the initial condition for 0 < x < 0.5, and R, for 
0.5 < x < 1. Figure 7.1 shows the results with 100 cells, h = 0.01, Xj = (j — 1/2 )h. The 
final time is t = 0.2 (after 57 time steps). The solid line represents the exact solution; 
the dots, the numerical one. Note that we do not need any boundary conditions because 
the final time is chosen so that the waves have not arrived at the boundaries. At the 
end of the calculation, there are a total of fourteen ‘nonsmooth’ cells: eight around 
the contact discontinuity, and six around the shock; the two discontinuities at the two 
ends of the expansion fan are in the ‘smooth’ region because the changes in density are 
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within the noise level by this time. The internal energy (per unit mass) e/p — u 2 / 2, 
which magnifies any error for this problem, is shown in Fig. 7.1(d). 

In the next two figures, the results for the density, velocity, and pressure are omitted 
since they are practically identical to those of Fig. 7.1. Figure 7.2 shows the internal 
energy of the same problem except for: (a) CFL number 0.9 (51 time steps), (b) CFL 
number 0.01 (4329 time steps), (c) k = 0 (57 time steps), and (d) k = 10 (also 57 time 
steps). As shown by Fig 7.2(a), in spite of the rather large CFL number, the solution 
is highly accurate thanks to the modification of the first ten time steps via expression 
(7.1). In Fig. 7.2(b), after a large number of time steps, the contact discontinuity still 
remains sharp due to the slope- steepening technique. Figure 7.2(c) shows that for Sod s 
problem, with CFL number 0.8, the final time is not large enough to smear the contact 
discontinuity because even without slope steepening (/c = 0), the solution is adequate. 
Comparing Fig. 7.2(d) to Fig. 7.1(d) we notice that once the contact discontinuity is 
resolved, increasing k does not yield a significant improvement. 

Figure 7.3 shows again the internal energy for the same problem except now we test 
the upwind fluxes: (a) central difference with artificial viscosity (5.13b) at ‘smooth’ re- 
gions, (b) one-intermediate-state (HLV) model everywhere, (c) primitive- variable split- 
ting (5.6) at ‘smooth’ regions, and (d) simple-wave upwind flux (5.75) everywhere. As 
shown in Fig. 7.3(a), the result by the central difference (5.13b) is as good as that 
via a more sophisticated model. However, for the Lax problem below, the solution by 
this model cannot match that via an upwind one. Figure 7.3(b) shows that the mod- 
ification (5.14) for the one- intermediate-state model is capable of resolving a contact 
discontinuity. Comparing Fig. 7.1(d) to Fig. 7.3(c), we notice that the solution via the 
primitive- variable splitting (at ‘smooth’ regions) is essentially identical to that via the 
one-intermediate-state model. Similarly, Figs. 7.1(d) and 7.3(d) show that at ‘smooth 
regions, a more sophisticated (simple- wave) model does not improve the result; in addi- 
tion, for a moving shock, the simple- wave model produces a result essentially identical 
to that via the shock model (flux-difference splitting). 

The next test, used by Lax (1954), is the Riemann problem 

(pl,ul,Pl) = (0.445,0.698,3.528), ( Pr,ur,pr ) = (0.5,0,0.571) (7.3) 

with the final time 0.16 (after 97 time steps). This problem is more difficult than that 


53 
















energy energy 





Figure 7.3. Sod’s problem with different upwind fluxes. 
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of Sod because the density field has to be ‘built-up’. The results are shown in Fig. 7.4. 

Figure 7.5 shows the density and velocity fields of the same problem except again 
(as in Fig. 7.3) we test the upwind fluxes: (a) ‘central difference’ (5.13b) at ‘smooth’ 
regions, and here note the bump at the tail of the expansion fan in the velocity field; 
(b) one-intermediate-model everywhere, and for this model, the bump is considerably 
smaller; (c) primitive- variable splitting at ‘smooth’ regions; and (d) simple-wave up- 
wind flux everywhere. 

The above two (and also the following) problems show that at ‘smooth’ interfaces, 
the one-intermediate-state model is as accurate as a more sophisticated model. It is 
also very economical. The flux-difference splitting with admissibility correction (or the 
simple-wave upwind flux) is needed only near a discontinuity. 

The third problem is designed to test the admissibility condition. It is the Riemann 
problem corresponding to a steady Mach 3 shock, except the left and right states are 
reversed: 

{pL,u L ,p L ) = (3.857,0.920,10.333), (pr,u r , Pr ) = (1,3.550, 1). (7.4) 

The solution for the final time t = 0.09 (after 59 time steps) is shown in Fig. 7.6. 

Note that for this problem, the shock for the numerical solution is about two mesh 
points behind that of the exact one. To explain this phenomenon, we plot the solutions 
in Fig. 7.7 for the density field at time t = 0.003 (after 4 time steps, and the expansion 
fan is well captured even here), t = 0.01 (9 time steps), t = 0.03 (23 time steps), and 
t = 0.04 (27 time steps). Because the numerical solution resolves the fan first, and 
then the contact discontinuity, by the time the shock is formed, it is already two cells 
behind. 

Also notice that if the flux-difference splitting without any correction is employed, 
then theoretically, the solution does not change with time because the fluxes are bal- 
anced. Eventually however, the solution is broken up by numerical error, but is no 
longer time-accurate. Next, a remark concerning the various admissibility conditions is 
in order. The correction for the flux- vector splitting in §5.5, that for the flux- difference 
splitting in §5.6, and the simple-wave upwind flux in §5.7 all produce essentially iden- 
tical results for this test. 

The fourth problem, due to Shu and Osher (1989), has several extrema in the 
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Figure 7.4. Lax’s problem. 
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smooth regions. It is used here to test the accuracy of the schemes. In the interval 
-5 < x < 5, a moving Mach 3 shock interacts with sine waves in density as described 


< - 4 , 

otherwise. 


(7.5) 


by the following initial conditions: 

_ / (3.857,2.629, 10.333) ifx 
\Pi u iP) | (1 + 0.2sin5x,0, 1) oth 

The final time is t = 1.8. Since the exact solution for this problem is unknown, the solid 
line is the solution with 1600 cells and k = 0. Figure 7.8 shows the solutions with 400 
and 800 cells, respectively (k = 5). Note that the solution for 800 cells has essentially 
converged (in the sense of mesh refinement). The solutions compare favorably with 
those of the third-order ENO scheme (Shu and Osher 1989). 

Our fifth test is the problem of two interacting blast waves (Woodward and Colella 


V = 


VI 

H 

VI 

1 is 

(Vl 

if x < 0.1, 

V M 

if 0.1 < x < 0.9, 

Iv* 

if 0.9 < x, 


(7.6) 


where 


Pl = PM = PR= u L = u M = u R = 0, PL = 10 3 , Pm = 10 2 , pr - 10 2 . (7.7) 


The boundaries at x = 0 (corresponding to the index 1 / 2 ) and x = 1 (index N + 1/2 
where N is the total number of cells) are solid walls with the following reflecting 
boundary conditions: at time level n (the superscript n is omitted below), we define 
auxiliary states V 0 , V-!, V_ 2 for the left boundary, and Vjv+i, Vn+ 2 , Vjv+3 for the 
right one by 

/?_;+! = pi, u-i+i = — ui, p-i+ 1 — pi, l — 1>2,3, (7.8a) 

Pn+i = PN-i+ii u n+i = — un-i+u Pn+i = PN-i+ii l — 1 , 2 , 3 . ( 7 . 86 ) 

At the final time t = 0 . 038 , the flow field has three contact discontinuities. The middle 
one, which is created by the two shocks interacting with each other, is very difficult to 
resolve. Even very compressive schemes such as Roe’s ‘Superbee ’ smear it (Roe and 
Pyke 1984 ). 

Notice that due to the small pressure in the middle region, the reconstruction step 
may result in negative values for pressure and density. This problem is solved simpl\ 
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Figure 7.8. Shu and Osher’s problem. 
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by replacing p by max(10 5 , p) (similarly for p) each time a new p or p is calculated 
(for ‘nonsmooth’ cells and interfaces). 

Figure 7.9 shows the solution with 1200 cells, k = 12, after 2411 time steps. In 
Figs. 7.9-10, this solution is represented by the continuous line since the exact solution 
is unknown. At the end of the calculation, there are nine ‘nonsmooth’ cells around 
each of the three contact discontinuities, and seven or eight around each of the two 
shocks. Note that again the solution compares favorably with those of higher-order 
schemes (Woodward and Colella 1984, Harten, Enquist, Osher, and Chakravarty 1987, 
Harten 1989). We also tested several values of k between 10 and 99; the results are 
practically identical for this problem. 

Concerning computing time, on an Iris Indigo R4000 workstation (16 Mflops, com- 
plied with option -02), this case takes 32 seconds. At ‘smooth’ regions, if the central 
difference with artificial viscosity (5.13b) is employed, there is roughly a 6 percent 
savings; if Roe’s flux-difference splitting (with no admissibility correction) is used, the 
cost is about 35 percent more. The three results are practically identical. 

Figure 7.10 show the density field of the solution with (a) 600 cells and k — 16; (b) 
400 cells and k = 20; (c) 400 cells and k = 7, where the middle contact discontinuity 
is slightly smeared, but the other two still remain sharp; and (d) 400 cells and k — 5, 
where the three contact discontinuities are smeared. 

The last example is (essentially) the quasi one-dimensional nozzle flow problem 
employed by Shubin, Stephen, and Glaz (1981): the domain is 2 < x < 10; the area, 
A{x) = 1.398-f .347tanh(.8x — 4); the supersonic inflow conditions are p = 1, u = 1.30, 
p — 1; the subsonic outflow, p = 1.4. The initial conditions are identical to the inflow 
conditions everywhere. Figure 7.11 shows the solution with 80 cells after 2000 iterations 
when machine accuracy is reached. Note that at the shock, the exact value for entropy 
calculated from the exact averages of density, momentum, and energy is plotted with 
the cross symbol. As a result, the entropy value at this point is meaningless. Also 
notice that the entropy field for this problem is difficult to resolve accurately; see, e.g., 
Hirsch (1990). 

Finally, we remark that for the case of a very slowly moving shock (Roberts 1990, 
Lin 1991, Quirk 1993), similar to the solution by the exact Riemann solver, all upwind 
fluxes presented here yield oscillatory results. We omit the details. As stated by Roe 
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(1989), “The perfect Riemann solver for generating numerical schemes probably does 
not yet exist.” 

8. Discussion and conclusions. The interpolation and constraints in this paper 
can be extended to the piecewise-parabolic, dual- variable (scheme V in Van Leer 1977), 
and dual- mesh methods (Nessyahu and Tadmor 1990, Sanders and Weiser 1992). The 
author hopes to report on some of the results in forthcoming papers. 

In summary, concerning the reconstruction step, a piecewise linear interpolation 
that preserves monotonicity and uniform second-order accuracy was introduced. The 
concept and coding of the monotonicity constraints were simplified by the use of the 
median function. Computational efficiency was enhanced by devising a criterion that 
detects the ‘smooth’ part of the data where the constraint is redundant. A slope- 
steepening technique, which can generally resolve a contact discontinuity in four cells, 
was presented. As for the upwind step, upwind fluxes were employed in a manner 
slightly different from those in the literature. They were rederived from a simple and 
unified point of view by approximating the Euler equations in conservation form using 
lin earization and diagonalization. At smooth regions, the one-intermediate-state model 
was employed for economic reasons. (The central difference with artificial viscosity 
(5.13b), which is even more economical, can also be used, but for the Riemann problem 
of Lax, the result is not as accurate.) Near a discontinuity, either the flux-difference 
splitting with admissibility condition or the simple-wave upwind flux was employed for 
robustness. The resulting scheme is efficient and accurate. 

For slow moving shocks, difficulties remain. Extensions to the multi-dimensional 
case, which may involve more theoretical work on constraints as well as upwind fluxes, 
are to be explored. 
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